Hybrid Equivalent Circuit/Finite Element/Boundary Element Modeling for Effective Analysis of an Acoustic Transducer Array with Flexible Surrounding Structures

: Transducer arrays are commonly analyzed by the finite element method (FEM) with high accuracy, but it is costly, particularly when having flexible surrounding structures. In this study, we developed an equivalent circuit (EC)-based model of an array of transducers with flexible surrounding structures for effective analysis. The impedance matrix was first constructed by coupling the electrical, mechanical impedance, and the acoustic radiation impedance obtained by the EC method and finite element-boundary element (FE-BE) coupling method. The transfer matrix of far-field pressure to the transducer response was then constructed by the FE-BE coupling method, and finally the sound pressure of the external acoustic field was obtained. To verify the accuracy, the results of the proposed method were compared with those of the conventional FEM. To evaluate the efficiency of the proposed method, the reduction in the degrees of freedom (DOFs) of the proposed method from the conventional FEM analysis was investigated. The simulation results of the proposed method are highly accurate and efficient. The proposed method is expected to be useful for conceptual design.


Introduction
Transducer arrays are usually applied for ultrasonic imaging [1][2][3][4][5] as they can detect objects in an acoustic medium by combining the electrical signals from each transducer.
Therefore, it is important to precisely predict the array performance while taking into account the multi-physical behavior in the transducer array.
FEM has been widely used for analyzing transducer array systems because of its high accuracy and the detailed information it provides for the design of transducer arrays. For example, Mestouri et al. [6] and Yamamoto et al. [7] used FEM to study the design of lowfrequency transducer arrays. Although FEM analysis is useful for final tuning, it is inefficient in the concept design stage because of its high computational complexity.
To help with this problem, simpler models of transducer array systems have been studied, such as EC models. The EC model of the transducer array can efficiently predict the multi-physical response of an array using a small matrix equation. There have been many studies on how to represent various types of transducers as EC models [8][9][10][11][12][13][14][15][16][17]. For example, Harrie extensively described the EC representation of various types of electromechanical transducers by using lumped parameter systems [8]. Later, he extended the theory and described the EC representation using distributed parameter systems [9].
In addition, the acoustic interaction between the acoustic field and transducer should be considered to apply the EC model to an array of transducers operating underwater, as in the case of sonar. The acoustic interaction effect is often characterized by the acoustic radiation impedance. For the calculation of acoustic radiation impedance, it is typical to use an analytical formulation or a numeric formulation. The analytical formulations have an advantage of computational efficiency. Many researchers have theoretically calculated the acoustic radiation impedance of pistons in an infinite rigid baffle [18][19][20][21][22][23][24][25][26][27][28][29][30]. Then, the underwater transducer array was tried to model using EC model with the acoustic radiation impedance of pistons in an infinite rigid baffle calculated by the theoretical formulation [10][11][12].
In many cases of actual sonar arrays, the acoustic radiation impedance usually cannot be obtained analytically as sonar arrays are finite and are sometimes not large compared to the wavelength in water [23]. To analyze the finite transducer arrays using the EC model, studies have used approaches based on Green's function to solve the Helmholtz integral equation [13][14][15][16]. For example, Audoly [13] applied this method to an array of circular pistons with a finite planar baffle to calculate the transmission characteristics of a sonar array. Yokoyama [14] applied this method to an array of rectangular pistons with a finite planar baffle to analyze the effects of acoustic interaction on the vibration velocity of transducers. Meynier et al. [15] developed a model based on this method to study the acoustic behavior of a capacitive micromachined ultrasonic transducer (CMUT) array.
However, the acoustic radiation impedance of an actual array having a complex shape is mostly coupled with the influence of surrounding structures. Actual sonars are affected by the acoustic effects of surrounding structures, such as baffles, enclosing domes, and nearby reflecting surfaces [16]. Exceptionally, Anthony [16] proposed an analysis method for array systems that reflects the influence of surrounding structures. He calculated the acoustic radiation impedance using only the Helmholtz integral equation formula, and so the surrounding structures were considered as rigid reflectors. Therefore, the influence of the vibration mode of the surrounding structure cannot be considered. In addition, this approach neglects the effects of sound waves that are transmitted through the acoustic window. Due to this limitation, studies on EC models so far have not been able to precisely predict the array performance while considering the flexibility of surrounding structures.
In this study, we developed an EC-based model of a transducer array system that can consider the effects of the acoustic interaction and the vibration of the flexible surrounding structures. With the proposed model, it is possible to precisely predict the response of the transducer array, including the crosstalk phenomenon. Crosstalk is a phenomenon in which the acoustic pressure generated by the projector is transferred to adjacent hydrophones through the surrounding structures, and the transferred pressure generates noise signals in the hydrophones. The performance of the transducer array is deteriorated due to this acoustic interaction [31].
The proposed analysis procedure can be roughly divided into two parts: One is the construction of the impedance matrix using the EC model and the FE-BE coupling model [32][33][34][35] and then using it to obtain the response of each transducer. To be more specific, the impedance matrix equation is derived based on the T-network EC model. The acoustic radiation impedance is calculated from the FE-BE coupling model. The FE-BE coupling model consists of the FE model of the surrounding structure, the BE model of the fluid between the array and the surrounding structure, and the BE model of the fluid outside the surrounding structure. The other part is the construction of an acoustic transfer matrix using the FE-BE coupling method and then the calculation of the sound pressure of the external acoustic field. With the proposed method, the array performance including the far-field directivity pattern and transmitting voltage response (TVR) can be estimated with less computational resources than the conventional FE analysis. Section 2 presents methods for constructing the impedance matrix equation and acoustic transfer matrix equation. In Section 3, a case of a cylindrical array with an enclosed dome is evaluated and compared to the result of multi-physics FEM to verify the accuracy and focus on a practical problem. We also evaluate the efficiency of the proposed method by comparing it with a conventional method based on Multiphysics FEM in terms of the reduction of DOFs. Conclusions are presented in Section 4.

Single Transducer
Tonpilz transducers are piezoelectric transducers mainly used for relatively low frequencies and high-power sound emission. Figure 1 shows a free body diagram of a tonpilz transducer. The transducer consists of its four basic components: the piston head mass, piezoelectric stack, tail mass, and stress rod. 1    The transducer's head mass is assumed to be rigid because it needs to provide uniform longitudinal motion, and its first flexural resonance should be significantly above the operating frequency band. Thus, the force 1,Head where j 1 = − , ω is the circular frequency, and Head M is the mass of the head.  1/ j where 0 C is the clamped capacitance, f C is the free capacitance, T ε is the permittivity coefficient, and N is the electromechanical turn ratio.
As shown in Figure 2b, the mechanical part of the EC model of a piezoelectric stack is represented using a T-network 12 based on a simple longitudinal bar model. The T-network is one of the ways to express a distributed parameter system as an EC model. The mechanical impedances of piezoelectric element 2 2 , c b Z Z can be written as where , c, A ρ , and L are the density, wave velocity, cross-sectional area, and longitudinal length of the part, respectively. k is the wave number. The subscripts PZT denote the piezoelectric stack.
The tail mass and stress rod have the same form of T-network EC model as the PZT stack, except that there is no electrical part. Thus, the following equations can be obtained in a similar way: The mechanical impedances of tail mass Z Z can be written as 3 Tail Tail Tail  Tail Tail   3 Tail Tail Tail  Tail Tail   4 s ts t s t s ts t 4 s ts t s t s ts t Here, the subscripts Tail, and st denote the parts of the transducer. The reaction force of the medium to the motion of the head surface r F is where S is the wet surface (the radiating surface) of the transducer head, r  is the position vector on the radiating surface, and p(r)  is the pressure produced by the transducer head. For a single transducer, r F can be expressed in terms of the self-radiation impedance r Z by dividing by the velocity of the head 1 u : By using the Equations (1), (2), (4), (7), and (10), the impedance matrix equation for a single tonpilz transducer can be expressed as Zu = F (11) where Z is the impedance matrix, u is a vector of velocities and current of the transducer, and F is a vector of the external forces and voltage of the transducer. They can be written as

Array of Transducers
For an array of transducers, depending on where the transducer is in the array, the reaction force of the medium, r F , has different characteristics. The force exerted on the i th transducer, , r i F , by the pressures from all the transducers is where n is the number of transducers that make up the array, i S is the wet surface of Here, , r ik Z is defined as the mutual radiation impedance between the i th and k th transducer, and 1,k u is the velocity of k th transducer head. Equation (15) shows that the coupling between transducers occurs in the acoustic region and can be expressed as mutual radiation impedance.
Generally, the transducers that make up an array are all the same transducers, so when constructing the impedance matrix of an array, all transducers have the same electrical and mechanical impedance. Thus, the impedance matrix equation of the i th transducer can be expressed as i i i Z u = F (16) where i Z is the impedance matrix of the i th transducer, i u is a vector of velocities and current of the i th transducer, and i F is a vector of the forces and voltage of the i th transducer. They can be written as where 1,i u , 2,i u , and 3,i u are the velocities, and i I is the current of the i th transducer.  (17), the impedance matrix equation for the array of transducers can be derived as where ik δ is the Kronecker delta, and ik Z  is an impedance matrix representing the relation between the i th transducer and the k th transducer.

Acoustic Radiation Impedance
In this study, an FE-BE coupling method has been developed to calculate the acoustic mutual-radiation impedance, , r ik Z , for an array of transducers with surrounding structures. This is done to reflect the influence of the vibration mode of the elastic surrounding structures and the influence of sound wave transmission and refraction by the surrounding structures on the impedance.
In this section, the formulations used are based on the coupling of the FE for the structure and the BE for the fluid. The FE/BE coupling method leads to reduction of the computational resources compared to when using the conventional FE model.
When an elastic structure of arbitrary shape is immersed in the fluid as shown in Figure 3a, the finite element modeling of the structure is given by the following discretized governing equation [34]: where M , C , and K are the mass, damping, and stiffness matrices of the elastic structure obtained by the FEM, respectively. s v is a velocity vector corresponding to all DOFs of the structure. s F is the harmonic load vector for the structure. G is a matrix that is used to find nodes on the surface (wet surface) where the elastic structure is in contact with the acoustic medium. Q is a transformation matrix that converts the outward unit normal force of the nodes on the wet surface to the force of the node in the fixed coordinate system (the global coordinate system). S is a diagonal matrix where the diagonal elements are the areas of each element on the wet surface of the elastic structure. p is the surface pressure acting on the radiating surface.
(a) (b) In this study, Equation (20) has been developed in the form necessary to obtain the acoustic radiation impedance of an array. When an elastic surrounding structure of arbitrary shape is placed around the radiating surface of the array as shown in Figure 3b, the finite element modeling of the surrounding structure is given by the following discretized governing equation: The subscript in denotes that the value is not on the wet surface. 0 1 2 s , s , s , and 3 s , indicating the sizes of each matrix, are the DOFs of the structure, the DOF of the fluid on the outer surface of the array, the DOFs of fluid 1 on the inner surface of the surrounding structure, and the DOFs of fluid 2 on the outer surface of the surrounding structure, respectively. Next, fluid 1 between the array and the surrounding structure and fluid 2 outside the surrounding structure satisfy the Helmholtz equation. The radiating surface is divided into boundary elements, and the fluid is modeled and formulated. Using the acoustic BEM, the relation between the surface pressure and the normal surface velocity of the radiating surface is given by the following equation: { } Assuming that there is no external force to excite the structure (

{ }
Now, where j S is the area of the j th boundary element.

Transfer Matrix for Field Point Pressure Calculation
The acoustic pressure field radiated from an array can be obtained by the acoustic transfer matrix and the velocity distribution over the transducer heads calculated by the impedance matrix equation in Equation (19). The transfer matrix, H , indicates the transfer function between the surface velocity of transducer heads, h u , and the pressure of field points, f P : The vector f P ( )

Analysis Model
The results of a general-purpose package (Multiphysics FEM software) were compared with the results obtained by the proposed method for a model with surrounding structures to confirm the validity of the proposed method.
The analysis model is a cylindrical array of 96 piston transducers arranged in 24 staves surrounded by an acoustic window, as shown in Figure 4. Figure 4a shows the geometry, and Figure 4b shows a FE-BE coupling model used in the proposed method to obtain the acoustic radiation impedance, , r ik Z , and the transfer matrix H . In Figure 4b, only the ellipsoid-shaped acoustic window is the structural finite element, the fluid 1 between the array and the acoustic window and the fluid 2 outside the acoustic window are the boundary elements. The input boundary condition of the FE-BE coupling analysis is shown in Equation (28). The acoustic radiation impedance, , r ik Z , can be obtained by postprocessing the analyzed surface sound pressure result as Equation (29). Although the transfer matrix H can be obtained using Equation (35), it is often difficult to access the system matrices (e.g., M , C , K , 1 A , 1 B , 2 A , and 2 B ) of the general-purpose package.
In that case, it can be obtained by setting the input boundary condition as Equation (28) and postprocessing the sound pressure result in the point of interest field. For example, the k th column of the transfer matrix H can be obtained by setting the input boundary condition as Equation (28). Figure 4c shows the conventional FE model. In the conventional FE model, an open boundary condition, such as perfectly matched layer (PML), is given to the surface of the sphere that is the external fluid domain. Moreover, each transducer is given a boundary condition that allows it to move in only piston mode without rocking mode because the transducer is designed so that the rocking mode does not occur.   Figure 5a shows the tonpilz transducer used in the array, and Figure 5b shows the transducer FE model used for the conventional analytical method. The material of the piezo-stack is PZT-4, and that of the tail and stress rod is AISI 4340 steel. The density of the head is 2700

Analysis Condition
This section describes how to determine the input vector i F of the impedance matrix equation. It is assumed that no external forces are applied to the array as the example analysis model is a projector array. The conditions are expressed by Equation (36).
Phase shifting the applied voltages is done to steer a projected beam. Normally, about 120° of the cylindrical surfaces is used for beam steering [17]. Thus, only 8 of the 24 staves are used in this model, which means that 32 of 96 transducers are operated, as shown in Figure 6.      Figures 8 and 9, respectively. The amplitudes of the transfer function between the applied voltage and the head velocities of transducers #7 and #9 are shown in Figures 10 and 11, respectively. Graphs are shown for elevation angles of the projected beam of 0° and 30°. Here, the analyzed frequency range is approximately

Analysis Result
The conventional FEM results and the results of the proposed method match well. The validity of the proposed method can be confirmed by the fact that the transfer functions of each transducer obtained by the two methods are very similar. In Figures 8-11, the analysis results are rather complex. This occurs because the sound diffraction, transmission, and vibration mode of the acoustic window have complex effects on the acoustic field produced by the sonar (the crosstalk phenomenon).

Array Performance
In this section, the beam pattern and TVR results calculated by the proposed method are compared with the conventional FEM results. The beam pattern and TVR are obtained by postprocessing the pressure at the field point of the far field. In this example, as the sound field at a distance of 100 m from the array has sufficiently developed into the far field, the field points are placed on a circle with a radius of 100 [m] on the y-z plane as shown in the Figure 12 where ( ) rms P r is the rms sound pressure measured at a distance r from the source along the acoustic axis. Here, r is multiplied to ( ) rms P r as TVR is normalized over distance r . Figure 12. Location of field points for calculating beam pattern.
The resulting TVR is shown in Figure 13. Results are shown for elevation angles of the projected beam of 0° and 30°. Typically, the operating frequency band of the array is designed as the frequency band near the maximum peak of the TVR. As shown in Figure 13, the maximum peak of TVR occurs between about The validity of the proposed method confirms that the array performances obtained by the two methods match. These analysis results for the exterior acoustic field are rather complex. In particular, the TVR results have several small peaks as shown in Figure 13. This occurs because the sound diffraction, transmission, and vibration mode of the acoustic window have a complex effect on the acoustic field produced by the transducer array (the crosstalk phenomenon).

Evaluation of Numerical Performance (in Terms of DOFs Reduction)
This section describes the numerical performance of the proposed method in terms of DOFs reduction. For numerical models such as FE model or FE-BE coupling model, DOFs is the number of unknown variables of the matrix equation for calculating the dynamic behavior of a physical system. The higher the DOFs, the greater the computational effort. Therefore, from the reduction rate of DOFs, we evaluated how high numerical performance the method proposed in this study has compared to the conventional full FEM.
We modeled cylindrical arrays as shown Figure 4a and configured four different sets of row and stave numbers, NC and NR to consider various sizes of problems: NC = 24 and NR = 4; NC = 36 and NR = 6; NC = 48 and NR = 8; and NC = 60 and NR = 10. We also considered four different analysis frequencies represented by ka: ka = 2, 3, 4, and 5. The following are basic schemes for generating the FE-BE coupling model for the proposed method and the FE model for the conventional method: 1. The size of one side of the acoustic element is set to about / 8 λ [36]. 2. The FE of the structure should be generated so that the mode shapes that are within 2 times the maximum analysis frequency can be well represented. 3. The shape of the acoustic window is assumed to be an oblate spheroid, and its size is set to increase in proportion to the array size. Here, the sizes of the half-length of the principal axes are assumed to be The DOFs for the two models (FE-BE coupling model and the conventional FE model) were obtained using the general-purpose software's DOF calculation function. The DOFs required for a reliable solution are summarized in Table 3. The DOFs of the proposed method are denoted by the broken line, and the solid line represent the cases of the transducer arrays corresponding to "NC = 24 and NR = 4", "NC = 36 and NR = 6", "NC = 48 and NR = 8", and "NC = 60 and NR = 10", respectively. As shown in Figure 16, higher frequency and a larger size of the array result in higher values of

Conclusions
This study developed a hybrid modeling method for an effective analysis of an acoustic transducer array with arbitrarily shaped flexible surrounding structures. More precisely, the proposed analysis procedure is divided into two parts: The first part is for the construction of the impedance matrix using an EC model and a FE-BE coupling model. We derived the electro-mechanical-acoustic coupling impedance matrix equation and established the procedure of calculation of the impedance matrices. Specifically, the mechanical and electrical impedances were derived from the LPM and T-network, a type of the EC model. The acoustic radiation impedances were calculated from the FE-BE coupling model, which represented the surrounding structure and acoustic coupling. The second part is dedicated to the construction of the acoustic transfer matrix using the FE-BE coupling method. Here, the acoustic pressure on each head surface can be calculated using the acoustic radiation impedance and the head velocities. The far-field acoustic pressure can be also evaluated by the acoustic transfer matrix equation and the head velocities.
The proposed method was applied to sonar with 96 tonpliz transducers and surrounding structures. The response of each transducer was calculated (such as the electrical admittance and the vibration velocity) along with the array performance (such as the farfield directivity pattern and the TVR). The results were compared with the results of FEM software. The results showed good agreement. To evaluate the numerical performance of the proposed method, the DOFs of the FE-BE models and the full FE model were investigated and compared. The DOF reduction rate of the proposed method to the conventional FEM is increased as the analysis frequency, array size, or size of the surrounding structures increased.
The main conclusion is that this proposed hybrid method can be used to effectively estimate the performance of the array while taking into account the surrounding structure and acoustic medium coupling. The computational effort is significantly reduced with the equivalent accuracy to the conventional FE analysis. Therefore, the proposed method is expected to be useful for conceptual design that require frequent design changes.