Abstract

A new approach is proposed for modelling medical ultrasonic transducers operating in air. The method is based on finite elements and the local interaction simulation approach. The latter leads to significant reductions of computational costs. Transmission and reception properties of the transducer are analysed using in-air reverberation patterns. The proposed approach can help to provide earlier detection of transducer faults and their identification, reducing the risk of misdiagnosis due to poor image quality.

1. Introduction

Ultrasonic transducers are widely used for medical applications. Reliable medical diagnosis requires good quality images from such transducers. However, it is well known that material degradation and possible faults to transducer’s elements often produce artifacts in ultrasonic images, as discussed in [1–4]. This can lead to false-positive and false-negative medical diagnosis [1, 5–8]. Therefore transducers need to be tested periodically. Professional systems (e.g., [1, 3, 4]) and tissue-equivalent phantoms are used in hospital practice to test transducers for image quality. Such tests usually performed in central hospital facilities are often not accurate, time consuming, and costly, as reviewed in [9, 10].

More recently a method based on in-air reverberation images, produced by a transducer operating in air, has been proposed [1, 8, 10]. In this test acoustic pulse is trapped in the transducers due to significant mismatch between acoustic impedance of lens and air. Reflections from various layers of the transducer produce characteristic reverberation lines in ultrasonic images, as illustrated in Figure 1. In-air reverberation patterns can be used to assess transmission and reflection properties of various transducer elements and their interfaces, which may be used to assess the image quality. Although the test associated with the entire procedure is relatively simple, interpretation of results and understanding of the method are still an issue. Firstly, current quality assurance testing, based on reverberation patterns, utilises only one simple parameter, that is, the distance to the last visible reverberation plane [5, 8, 11]. Clearly, this approach is not sufficient for accurate image quality assessment. Secondly, even if more parameters were used interpretation of results would not be an easy task due to the fact that there is very little knowledge of correlation between features in reverberation patterns and possible transducer faults. It is clear that efficient modelling and numerical simulations could provide valuable information on this correlation helping to provide better understanding of in-air reverberation patterns. As a result, earlier detection of transducer faults could be possible, preventing possible medical misdiagnoses.

Various approaches have been used for modelling of ultrasonic transducers, including analytic, semianalytic, and numerical techniques [12]. It appears that mechanical and equivalent electrical circuit coupled oscillators (e.g., electrical circuit oscillators [13–16]) and methods based on Finite Element (FE) analysis [17–25] are the most widely used approaches for medical transducer modeling. The former approach is relatively simple and efficient but provides only basic information on transducer properties. FE analysis is well suited for arbitrary geometries and, with the arrival of multi-physics software packages (e.g., PZflex and COMSOL), convenient for modeling of complex anisotropic and functional materials. This is essentially required for modelling of ultrasonic transducers and imaging arrays, as explained in [24]. Finite element analysis can provide detailed information on transducer performance but is computationally expensive. Previous work in this area has demonstrated that numerical simulations of ultrasonic pulse propagating within layers of a medical ultrasonic transducer exposed to air produce in-air reverberation pattern images that are similar to those obtained in hospital practice [8, 26]. These exploratory investigations have also showed that simulated in-air reverberation patterns can uniquely reveal faults in different elements of transducer. The major drawbacks of the work presented in [8] relates to computational efficiency. Even current two-dimensional (2D) numerical simulations of FE-based [26] wave propagation in ultrasonic medical transducers have led to very large and computationally expensive models.

The paper proposes a new approach that can be used for modelling of medical ultrasonic transducers operating in air. This approach is based on a hybrid methodology. The crystal is modelled using the classical FE analysis based on the constitutive electromechanical equations for piezoelectric materials. Wave propagation in backing, matching, and lens layers of the transducer is modelled using the local interaction simulation approach (LISA) that has been originally developed for wave propagation modeling in complex media with sharp interfaces [27–29]. There are two important motivation factors that have driven the research work presented in the current investigations. Firstly, the LISA has proven performance superiority over classical FE analysis in numerical simulations involving large impedance mismatches [30]. Secondly, recent developments in the LISA implementation, reported in [30, 31], have also proven computational superiority over wave propagation FE-based analysis. This new implementation of the LISA relates to a parallel algorithm for iterative wave propagation equations based on a new Computer Unified Device Architecture (CUDA) computation strategy, available in low-cost graphical cards [30].

The paper starts with a brief introduction to medical ultrasonic transducers. This introduction is given in Section 2. The theoretical background related to the proposed methodology is given in Section 3 for the sake of completeness. This includes the electromechanical constitutive equations and the LISA algorithm formulation, needed for wave propagation in the crystal and the remaining transducer layers, respectively. The proposed hybrid methodology is implemented in Section 4. Numerical implementation of wave propagation in one- and fifteen-finger medical ultrasonic transducers is presented. Numerical simulation results are given in Section 5. These results are compared with the classical FE-based approach. Finally, the paper is concluded in Section 6.

2. Medical Ultrasonic Transducer

A typical medical ultrasonic transducer, illustrated in Figure 2, is an array transducer that consists of four major elements. These are crystal (ceramic), backing, matching, and lens layers. The matching layer guarantees that acoustic impedance of the transducer is close to the impedance of the tissue. The function of the backing layer is to damp vibration of the ceramic in the direction opposite to the wave propagation direction. This is required to reduce undesired reflections and interferences in ultrasonic images.

The crystal (piezoelectric material) is a core element of the transducer. The inverse piezoelectric effect is used in the ceramic to generate ultrasonic waves. These waves travel through the matching and lens layers and then interact with human tissues. Wave reflections from various layers and tissues travel back to the piezoelectric material. The direct piezoelectric effect is then used to obtain voltage signals that are used to produce ultrasonic images. Amplitude peaks in these images indicate reflections from various layers and features. When speed of sound is known and peak arrival times are established, the relevant distances can be estimated to produce ultrasonic images for medical diagnosis purposes.

3. Hybrid Modelling for Transducers Operating in Air

This section presents the hybrid approach to modelling the medical ultrasonic transducer. Firstly, the constitutive equations for the piezoelectric crystal are presented for the sake of completeness. Then the local interaction simulation approach is briefly described.

3.1. Constitutive Equations for the Crystal Element (Ceramic)

The inverse and direct piezoelectric effects are used in the ceramic element to generate and sense ultrasonic waves, respectively. The inverse piezoelectric effect is the ability of material to deform under applied electrical field, whereas the direct piezoelectric effect is the ability of material to generate electrical charge in proportion of applied force. Both effects can be described using the constitutive equations that couple the mechanical and electrical fields [32]; that is, where is the electric displacement, represents the strain field, is the stress-based piezoelectric matrix, is permittivity, and is the electric field. The superscript in the above equations indicates the transpose and the superscripts and represent constant stress and electric fields, respectively. Therefore is the permittivity for constant mechanical stress and is the compliance for a constant electrical field with this respect. It is important to note that the Voigt notation was used in (1), so the compliance , normally a rank-four tensor, is a 6 × 6 matrix.

Two major simplifications to (1) can be made according to the dominating piezoelectric effect of the transducer, that is, transversal () and longitudinal () effects. For the piezoelectric transversal case, due to the dominating coefficient, it may be assumed that the only component of induced stress is perpendicular to the electric fieldFor the piezoelectric longitudinal effect the induced stress is parallel to the electric field since the dominating piezoelectric coefficient is . Then the stress field may be expressed asIn order to solve the entire problem, the electric balance, that is, and appropriate boundary conditions are needed additionally. The strain and electric matrices can be expressed in terms of mechanical displacements and voltage asEquations (1) to (5) are solved using the FE to find the displacement field of the piezoelectric crystal under applied voltage and the voltage under applied displacements for the actuation and sensing, respectively.

Numerical simulations with finite elements can be used not only to model the piezoelectric effects but also to investigate wave propagation problems in medical ultrasonic transducers, as illustrated in [8, 24, 26]. However, this approach is not effective when accuracy and computational effort are considered. The application of the LISA can offer a solution to the problem. Previous research works [27–31] show that the LISA is very well suited to model wave propagation in complex media with large impedance mismatches. Also, recent investigations, reported in [30], demonstrated that parallel architecture of the LISA implemented on graphical cards can significantly reduce the computational effort. These two important developments point to the application of the LISA for wave propagation modeling in medical ultrasonic transducers.

In the LISA approach the layers comprising an ultrasonic transducer are discretised into grid of cells. For each layer different material properties are assigned. Material properties are assumed constant within each cell but may differ between cells. The entire algorithm can be derived from the elastodynamic wave equation for homogenous media. This equation can be given as where and are Lamé constants, is the material density, is the vector of particle displacements, and a comma preceding a subscript denotes differentiation with respect to that variable. For a given material the Lamé constants can be expressed aswhere is the Young modulus and is Poisson’s ratio.

For the 2D wave propagation case investigated in this paper, the antiplane shear waves can be ignored and (6) can be simplified asand further rewritten in a matrix form aswhere the relevant matrices, containing material data, can be expressed asWhen the structure is discretised into cells, as illustrated in Figure 3, the junction of the four cells defines the nodal point P. The sharp interface model (SIM) condition is applied to derive the LISA interaction equation; it means each cell is treated as discontinuous and stresses are matched at interface grid points. The second time derivatives across the four cells are required to converge towards a common value at this point:This ensures that if the cell displacements are continuous at P for the two initial times, that is, and , they will remain continuous for all later times. A finite difference scheme is used to calculate spatial derivatives from (9) in the four surrounding cells to P. This leads to four equations by using (9) and (11):The subscripts denote cell number to which the data belongs, and an index after a comma denotes differentiation with respect to this quantity.

Figure 3: Discretisation scheme used for wave propagation modeling in the 2D LISA.

Since only four equations of type (9) are available, as shown in (12), additional four conditions must be found to solve for the displacement field of the structure. For this purpose the stress continuity across the cells’ interfaces is established. This leads to the set of eight equations in eight unknown quantities.

The components of the stress tensor in 2D for orthotropic materials [27] are given bywhereAs is pointed out before the material matrices , , , , and are constant inside each cell. Imposing continuity of the stress tensor at the point P and using further finite difference formulae give four additional equations in the unknown quantities , and :thus allowing these unknown spatial first derivatives to be solved. For more details the reader is referred to [27]. As a consequence of the SIM approach, the density , shear modulus μ, and were redefined at the interface as the averaged values taken from the neighboring cells. Therefore, this model can be used for both homogeneous and heterogeneous media. For simplicity some additional definitions are made: The particle displacements in the and directions are denoted as and, respectively.

Recalling that , after a certain amount of algebra, these and displacements at the point P can be calculated aswhere , . Equations (17) are the principal displacement equations of LISA in two dimensions. A more detailed derivation can be found in [27, 30].

It should be also pointed out that the ultrasonic medical transducer is a multilayered structure. The wave propagating through a single finger (one element) is partially reflected from the interfaces and fully reflected from the lens surface and acquired by the piezoelectric crystal forming the air reverberation pattern. Since this pattern is used to detect and identify damage of the MUT, particular attention should be paid to the wave interaction with materials’ interfaces. It was already mentioned that the LISA iteration equation involves stress matching across cells. This leads to the so-called sharp interface model (SIM) that is capable of solving wave propagation problems for sharp impedance mismatches very accurately. This allows for efficient damage and traction-free surface modeling that is used in this study.

4. Numerical Implementation

A 2D model, presented in Figure 4, of the medical ultrasonic transducer operating in air was used in numerical simulations. The four-layer structure (backing, ceramic, matching, and lens) is called a “finger” in this paper for simplicity; the study involved two simple configurations, that is, one- and fifteen-finger transducers.

Figure 4: Multi-finger model of ultrasonic transducer.

Numerical implementation described in this section relates to two major tasks: (1) ultrasonic wave generation by the piezoelectric crystal and (2) wave propagation in the remaining layers of the medical ultrasonic transducer. This problem is solved using a hybrid approach. Wave generation by the crystal (task 1) is simulated using finite element modelling and wave propagation in a multilayered structure (task 2) is simulated using the LISA-based modelling, as described below.

4.1. Ceramic Modelling—Ultrasonic Wave Generation

The constitutive electromechanical equations, described in Section 3.1, can be solved through approximations and discretisation of the transducer geometry. Finite element modeling, based on ANSYS, was applied. The coupled-field quadrilateral PLANE13 solid elements were used to model the piezoelectric layer of the transducer. These elements are defined by four nodes and three degrees of freedom per node (i.e., two displacements and voltage).

One 0.25 × 12 mm crystal finger (element) was modelled using 1200 finite elements and 1464 nodes. Large numbers of element were selected to avoid nonphysical frequencies and direction dependent dispersive character, following recommendations given in [8, 26]. The PZT-5A piezoceramic was selected for the material of the ceramic. Electromechanical properties of this material are given in Table 1. The thickness of the ceramic layer was defined by the required operation frequency of the transducer. To maximize the efficiency piezoelectric crystals are usually made with a thickness of around 1/2 of the demanded wavelength.

A pulse of one sine cycle of 4 MHz frequency and 300 V amplitude was used as the input signal to the piezoelectric crystal. When receiving signals were gathered, the upper electrode was grounded, and the voltage signal was collected from the lower electrode.

Once the geometry was discretized, material properties were selected, and excitation was defined, the solution was obtained for nodes of the elements. The required solutions for arbitrary positions were obtained using a linear combination of polynomial interpolation functions, as described in [33–36]. The entire analysis was performed in the time domain using a step-by-step time explicit integration scheme to obtain displacements. The time-step was always selected to achieve the Courant-Friedrichs-Lewy (CFL) stability condition. The transient analysis was performed from 0 to 5 μs with the interval of s and from 5 to 10 μs with the interval of s. The calculation time required to model one-finger transducer using FE was equal to approximately 15 minutes.

4.2. Wave Propagation Modelling

Since a layer of backing material is placed on all but one side of the transducer ceramic, the elastic waves generated by the transducer propagated in that one direction only. The backing layers damp out waves propagating in other directions. Therefore the displacement field of the ceramic was not distorted by waves reflected from boundaries.

The heights of the backing, matching, and lens layers were equal to 7.5, 2 and 2 mm, respectively. The width of all these elements was the same as the width of the ceramic. Discretisation in space and time was applied for wave propagation modelling. The transducer was discretised into 0.05 × 0.05 mm elements in the LISA model. The total number of elements used to model wave propagation in the fifteen-finger transducer was equal to 25 758. The time step used was equal to 0.005 μs to maintain the stability of the explicit time integration procedure. Once the model was discretised, material properties for transducer elements were selected. Backing layers are usually selected from a material that meets two conditions: (1) the ability to attenuate sound waves should be relatively high; (2) the acoustic impedance of breaking layer should be similar to the acoustic impedance of the ceramic layer for avoiding wave reflections from the ceramic-backing layer boundary caused by the acoustic impedance mismatch. Therefore the Tungsten-Epoxy material was used to model the backing layer. The matching layer was modelled using Araldite properties. The lens layer was modelled using properties of Perspex. The relevant material properties of the matching, backing, and lens layers are layers given in Table 2 as reported in [8, 26]. It is important to note that the piezoelectric crystals were not modeled explicitly with the LISA. The FE displacement field in the upper layer of ceramic was passed to the matching layer to provide excitation for wave propagation. The calculation time required to model one-finger transducer using LISA was equal to around 1 minute.

Table 2: Material properties used for transducer modelling.

The LISA wave propagation equations, described in Section 3, form a series of iterative equations of the formwhere denotes vector of field variables and is the time. Since the further time step is being calculated as a combination of the same quantities taken in already calculated time steps, the algorithms can be parallelised very efficiently, as described in [30, 31]. Therefore the LISA-based wave propagation modelling was implemented using the CUDA and run on Graphical Processing Units (GPUs). This low-cost solution led to significant reductions of computational costs. The total computation time to model the fifteen-finger ultrasonic transducer was reduced from approximately 4 hours [26] for the FE model to approximately 14 minutes for the combined FE and LISA model. For more detailed comparison of LISA model with traditional method like FE, the reader is referred to [37].

5. Numerical Simulation Results

Once the FE- and LISA-based model of the medical ultrasonic transducer operating in air was implemented, numerical simulations were performed to analyse transducer responses. These responses were captured at the bottom of the matching layer. The results for one- and fifteen-finger transducer models are presented in this section. These results (corresponding to the LISA solution) are compared with numerical simulation results based on FE analysis, where both actions, that is, ultrasonic wave generation and wave propagation, were modelled using FE.

5.1. One-Finger Transducer Model

The LISA and FE displacement responses for the one-finger (one element) transducer are presented in Figure 5. The results display three major wave packages indicated by arrows. These wave packages are reflections that can be identified knowing the relevant material properties and wave speeds. The first dominant wave package, observed after 0.22 μs, is the excitation signal. The second wave package, after 1.52 μs, corresponds to the reflection between the matching and lens layers. Finally, the reflection from the lens-air boundary can be observed as the third eave package after 4.1 μs. When the relevant amplitudes and arrival times are compared, the LISA-based results match the FE results relatively well. Nevertheless some differences can be observed, particularly when the 2nd and 3rd wave packages are compared. The LISA wave propagation response displays slightly larger values of arrival time (for 2nd, 6.5% and for 3rd, 2.5%) and smaller values of amplitude (for 2nd, 42.9% and for 3rd, 35.7%) for these wave packages.

After the relatively good agreement was obtained for the LISA- and FE-based numerical simulations, structural damage was introduced to the models. Damage to the lens layer was simulated by 50% stiffness reduction. Figure 6 gives the LISA and FE displacement responses for the one-finger damaged transducer. Although both responses display the same wave packages, as expected, some difference can be observed when the amplitudes and arrival times of the second and third wave packages are compared. The effect of damaged lens can be analysed when the displacement from the healthy and faulty transducers is compared in Figure 7. Here, the relevant displacements were obtained using the LISA-based numerical model of wave propagation. The results show that the damaged lens increases significantly the amplitude of the second wave package due to the impedance mismatch caused by stiffness reduction.

5.2. Fifteen-Finger Transducer Model

This section presents the results for the fifteen-finger transducer model. Figure 8 gives the displacement obtained from the LISA and FE models for the 2nd, 5th, 11th, and 15th finger of the transducer. The results display very good agreement between the FE- and LISA-based models of wave propagation. Slight differences can be accredited to the different formalism of the methods applied.

Figure 8: Medical transducer displacement responses for the fifteen-finger transducer. The responses are given for the selected transducer fingers, as indicated.

Numerical simulations were also performed for the damaged fifteen-finger transducer. The damage was introduced to the sixth finger and simulated by the 50% reduction of stiffness to the matching layer. The results for the LISA model are in very good agreement with the FE model, as shown in Figure 9. The effect of damage in the LISA model can be observed when the relevant displacements for the healthy and faulty transducer are compared in Figure 10. The wave reflection from the damage of matching layer can be observed between 0.22 and 1.52 μs. The former corresponds to the time of reflection from the boundary between the crystal and the matching layer corresponds to the reflection from the boundary between the matching layer and the lens.

Figure 9: Medical transducer displacement responses for the fifteen-finger faulty transducer. The response is taken form the sixth finger of the transducer.

Figure 10: Medical transducer LISA-based displacement responses for the fifteen-finger healthy and faulty transducer. Transducer responses are given for the 6th finger. Major differences indicated by the window.

The piezoceramic crystal displacement responses were transformed to the in-air reverberation image patterns, as described in [8, 10]. The reverberation patterns, obtained from the LISA model, for the healthy and faulty transducer are given in Figure 11. Here, the brightness in the image is proportional to the amplitude of the displacement generated by the piezoelectric crystal when a reversal pulse excites the ceramic. When the stiffness of the matching layer is reduced the displacement from the ceramic is affected, as expected. This can be observed also in the reverberation pattern in Figure 11(b). The results also show that the displacements of the two neighbouring fingers are also affected. The effect of damage, observed in the in-air reverberation pattern, is indicated in Figure 11(b).

6. Conclusions

A new, hybrid method was proposed for modelling medical ultrasonic transducers operating in air. This method is based on finite elements and the local interaction simulation approach (LISA). The former is used for ceramic modelling, whereas the latter is applied for wave propagation analysis. Wave propagation analysis was implemented on GPUs in order to reduce the computational effort. Numerical simulations were performed for the one- and fifteen-finger transducer and compared with the classical FE results.

This exploratory study shows that LISA produces very similar results to FE analysis. Small differences in amplitude and arrival times of reflected wave packages, observed for the one-finger transducer, disappeared when the fifteen finger transducer was analysed. The results also show that the LISA displays the same damage-related features in models of faulty transducers as the FE model.

Significant reduction of computational time, from days to minutes, is the major advantage of the proposed method. This is particularly attractive for future developments in which real medical transducers (e.g., based on 128 fingers) should be modelled in 3D configuration.

In summary, the work presented in the paper demonstrates that in-air reverberation images can be considered for fault detection in medical transducers. Efficient numerical simulations, based on LISA, can significantly ease the modelling task that is needed to implement the method in practice. Further research work is required to correlate features in in-air reverberation patterns with real transducer faults.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The work presented in this paper was supported by the Foundation for Polish Science under the research WELCOME Project number 2010-3/2 (Innovative Economy, National Cohesion Programme, and EU).