3-D Modelings of an Ultrasonic Phased Array Transducer and Its Radiation Properties in Solid

Ultrasonic phased array technology has become popular in the field of industrial nondestructive testing (NDT). The technology is notable for its capability to provide images of the inside of a target in real time. Since the image quality results from rapid steering and focusing of the radiation beam from an array transducer, it is essential to have well understood characteristics of the radiated beam. Even though fundamental concepts of the array transducer in medical fields were introduced over 30 years ago (Macovski, 1979), the beam modeling and optimization of the array transducer for the NDT have only been investigated intensively in the past decade (Azar et al., 2000; Song & Kim, 2002). The characteristics of the radiated beam from the ultrasonic phased array transducer vary according to the transducer parameters such as frequency, aperture size, number of elements, pitch(inter-element spacing), element width, layout dimensions and so forth. If these parameters are not chosen properly, spurious grating lobes or side lobes with high amplitude will exist in the radiation beam field. Consequently, the image quality will be deteriorated significantly.


Introduction
Ultrasonic phased array technology has become popular in the field of industrial nondestructive testing (NDT). The technology is notable for its capability to provide images of the inside of a target in real time. Since the image quality results from rapid steering and focusing of the radiation beam from an array transducer, it is essential to have well understood characteristics of the radiated beam. Even though fundamental concepts of the array transducer in medical fields were introduced over 30 years ago (Macovski, 1979), the beam modeling and optimization of the array transducer for the NDT have only been investigated intensively in the past decade (Azar et al., 2000;Song & Kim, 2002). The characteristics of the radiated beam from the ultrasonic phased array transducer vary according to the transducer parameters such as frequency, aperture size, number of elements, pitch(inter-element spacing), element width, layout dimensions and so forth. If these parameters are not chosen properly, spurious grating lobes or side lobes with high amplitude will exist in the radiation beam field. Consequently, the image quality will be deteriorated significantly.
In this chapter, we first show principles of electronic scanning with phased array transducers for the NDT and effective settings of the delay time for the beam steering and focusing. And then mathematical models of linear and matrix array transducers and simulation tools to predict ultrasonic beams from the transducers are outlined. We explain three-dimensional (3-D) numerical simulation tools in both frequency and time domains. This chapter is arranged as follows: •I n S e c t i o n 2, we show various ultrasonic array transducers for the NDT and electronic control of the array transducer. First the appropriate delay time setting for beam steering and focusing is described, and then causes and prevention of the grating lobe are explained.
•I n S e c t i o n 4 the finite integration technique (FIT) (Fellinger et al., 1995) is introduced as a time domain calculation tool. The FIT represents very stable and straightforward schemes to investigate the wave propagation in complex fields such as inhomogeneous and anisotropic media. For example, a dissimilar welding part includes metal grains of different size and orientation, therein the FIT is of assistance to predict the ultrasonic propagation direction as well as scattering and attenuation in the welding part. Here a transient wave simulation of the array ultrasonic testing (UT) for a welded T-joint is demonstrated.
•I n S e c t i o n 5, we summarize our research and discuss prospects for array UT.

Array transducer for NDT
Ultrasonic phased array techniques have been applied to the NDT. In the phased array techniques, angled beams with various focal lengths are generated by the array transducers composed of small ultrasonic transducer elements. The type of beams is electrically controlled through delay time control for excitation of transmitted waves and signal processing of received waves.
The array transducers have two advantages for NDT in comparison with the conventional single or dual element transducers. Firstly, the various angled beams can cover wider range of an inspection target. Secondly, NDT results can be visualized immediately as a cross section images (B-scan) or 3-D images.

Principle of electric scanning
General configurations of the array transducers are classified into one-dimensional (1-D) arrays and two-dimensional (2-D) arrays shown in Fig.1. Linear arrays ( Fig.1(a)) and annular arrays ( Fig.1(c)) are typical configurations of the 1-D array transducers. The linear arrays provide beam-steering and focusing in a plane. The Annular arrays adjust focal length along an axis. In recent years, the 2-D arrays have been applied to the NDT field for maximization of the advantages of the ultrasonic phased array technology (Drinkwater & Wilcox, 2006). Matrix arrays ( Fig.1(b)) and segmented annular arrays ( Fig.1(d)) can steer beams in 3-D volume. Ultrasonic beams are electrically scanned by delay-time control for each transducer element composed of array transducers as shown in Fig.2. Patterns of the delay time for transmitting and receiving (hereinafter referred to as "delay law") are stored on the delay-time controller which shifts excitation pulses in transmitting and signals of reflection waves in receiving.  The main parameters of the linear and matrix arrays are the aperture, the element pitch, and the element number, shown in Fig.3. The apertures A are expressed as a function of the element pitch p, the element number N and the element gap g in Equation (1). For the matrix array transducers, subscripts of parameters correspond to the x and y axes. Electrical scanning patterns widely-used in the NDT field are sectorial and linear scanning as schematized in Fig.4. B-scan results for electrically discharged notch of the sector and the linear scanning are displayed also in Fig.4. The through-wall depth of the notch is measured as the difference of the z-axis between echoes from the root and the tip of the notch. NDT data in the sectorial scanning is displayed as a sector-form B-scan image (cross section) by changing propagation directions. In the sectorial scanning, imaging area can be wider than an aperture of an array transducer. Element number is generally from 15 to 63. The linear scanning provides a parallelogram B-scan image by switching positions of active elements. In the linear scanning, imaging area can be wider by a large number of transducer elements, e.g. 63 to 255. Substituting a motion axis of mechanical scanning for the linear scanning can shorten scanning time of an automated UT.

Delay setting for electrical scanning
The coordinate systems for the linear array transducer and the matrix transducer are defined as Fig.5. The delay laws imposed on each transducer element positioned at x =( x n , y m ) for an angled beam with focal point F in a material characterized by velocity c is expressed by the difference of the propagation time: where x 0 is the center position of the array transducer. In this section, the center position is set to the origin. The focal point F is expressed by the incident angle (zenith angle) θ and rotational angle (azimuth angles) φ : F =( F 1 , F 2 , F 3 )=( R sin θ cos φ, R sin θ sin φ, R cos θ). Therefore the delay law ∆τ nm is obtained as In the case of the linear array transducers, Equation (3) is simplified to In the finite limitation of the focal length R, Equations (3) and (4) are simplified to the following expressions, respectively:

Effects of beam focusing and steering
Characteristics of the focusing beams , such as a beam width, peak amplitude and focus length, depend on both parameters of the array transducers (shown in Fig.3) and delay laws. As examples, delay laws of normal and angled focused beams for the linear array transducers are plotted in Fig.6. The horizontal axis is position of the transducer element and the vertical axis is delay time. Negative values of the delay time correspond to transmitting ultrasonic beam before the relative origin of the time. Therefore the delay time in an actual ultrasonic phased array equipment is shifted to nonnegative value by adding constant values to the theoretical delay time in Equation (3).
Beam profiles for the delay laws in Fig.6 are calculated based on a Rayleigh-Sommerfeld model (Kono et al., 2010), as shown in Fig.7. The calculation conditions are as follows: the array configuration is the linear array transducer with 16 elements; element pitches, a half of a wavelength λ of 2.95mm; an element width, 8λ; center frequency, 2.0MHz; focal lengths, 20-80mm and incident angles, 0-45 • . Peak levels of the amplitude have a propensity to decrease with increased focal length from the results of the normal beams in Fig.7(a)-(c). For the angled beam, peak levels decrease with higher incident angles from the results in Fig.7(b),(e) and (f).

Generation of grating lobes
Grating lobes can be generated under the condition that difference of the delay time between the adjacent transducer elements is a multiple of a period T. For the matrix array transducers, Equation (5) is separated into x-dependent and y-dependent terms: where the position of the transducer element is set to element pitches: x =(p x , p y ).Theprime mark indicates the direction of grating lobes, m and n are integer indexes for the grating lobe with phase difference 2(i + j)π.
Therefore, the propagation direction of grating lobes for the matrix array transducers are calculated from the delay time differences in Equations (8) and (9) (Kono & Mizota, 2011).
Similarly, the generation conditions of grating lobes for the linear array transducers are derived from Equation (11) substituted p for x n : where p is an element pitch and i is an integer index for the grating lobe with phase difference 2iπ. The incident angle of grating lobes for the linear array transducer is transformed from Equation (12) and the relation between cycle T and wavelength λ : T = λ/c :  Equation (13) is consistent with Equation (11) substituted zero for the rotation angles of the main lobe and grating lobes. Grating lobes are generated on the condition that the absolute value of the left side in Equation (13) is less than or equal to one; therefore the grating lobes can be suppressed when an element pitch is less than one-half wave length.
Beam profiles for angled beams generated by two kinds of linear array transducers with different element pitch are compared in Fig.8. The position of the focal point F is the same each other. No grating lobe is observed for the linear array transducer with element pitch equal to one-half wavelength. A grating lobe is generated in the case of the array transducer with one-wavelength element pitch. The incident angle of the grating lobe is obtained as −17.0 • from Equation (13) for i=1 and θ=45 • .

Modeling of phased array transducer in frequency domain
Fundamentally, ultrasonic beam can be modeled by superposing spherical waves according to the Huygens' principle. The Rayleigh-Sommerfeld integral (RSI) model is also based on the Huygens' principle and well known as an expression of wave field in liquid from an immersion piston transducer. On the other hand, the wave field in solid from a contact piston transducer can be expressed with the Vezzetti's model (Vezzetti, 1985). This is a little bit complicated because the Vezzetti's model remains the integral form of the fundamental solution in elastic half space. Schmerr (Schmerr, 1998) has introduced an explicit expression of the Vezzetti's model using the stationary phase method. Here we show numerical modelings of the ultrasonic phased array based on the Schmerr's expression.
Here we assume the harmonic wave of exp(−iωt) time dependency. First of all, consider a single planar piston transducer which generates a P-wave in solid. To model the radiation field from the transducer, we use the geometry shown in Fig.9, where the solid region is the half-space x 3 ≥ 0. On the x 1 -x 2 plane, we assume that velocity in the x 3 direction is zero everywhere except for a finite region S. Here the displacement in solid due to the P-wave is written as (Schmerr, 1998) u where, r = |x − y|, c p and k p (= ω/c p ) are the velocity and the wave number of the P-wave, respectively. Also ρ is the density and P 0 is a constant known pressure in S. In Equation (14) d p is the polarization vector: and K p (θ) is the directivity function: where c s is S-wave velocity, κ = c p /c s and θ is the angle between the vector (x − y)a n dt h e normal of the transducer surface.
In the next section, we show two methods to obtain the solution u in Equation (14). One is based on a numerical integral method and the other is an approximation method using superposition of the Gaussian beam. The former method needs some numerical calculation but is effective even if the transducer element has a complicated shape. The latter case is a fast method, but only applicable to circular and rectangular element of the transducer.

Rayleigh-Sommerfeld numerical integral (RSNI)
It is not easy to solve the integral in Equation (14) analytically except for simple element shape such as a rectangle or a circle. Therefore let us consider a numerical integral over S in Equation (14). For convenience we use the local coordinate system (ξ, η) on the transducer surface.
Using the four-node two-dimensional element, we can interpolate the coordinate y as where In Equation (17), y α is the coordinate of the four element nodes. Substituting Equation (17) into (14), we can rewrite the integral form as follows.
where |Y| is the Jacobian operator:  where w is the integration weight and I is the number of the integration points.
In the above procedure, the wave field from a piston transducer, namely the wave field from an element of the array transducer, can be calculated. When the array transducer has elements N in the x 1 direction and M in the x 2 direction, we can obtain total wave field by appropriately adding up N × M wave components. Using the time delay ∆τ nm for the element positioned at n-th in the x 1 direction and m-thinthex 2 direction, the total wave field can be expressed as Equation (22) can be applied to not only a flat array element but also a curved one.

Multi-Gaussian beam (MGB) model
The multi-Gaussian beam (MGB) model (Schmerr, 2000) is adopted for simulation of beam fields radiated from circular and rectangular transducers. However, the MGB model is based on the assumption of the paraxial approximation, and will lose its accuracy at the high-refracted angle, which is a condition that paraxial approximation is not satisfied (Park et al., 2006). Such problem often happens in simulating the steering beam of phased array transducers with the conventional MGB model. To address such a problem, the MGB without the nonparaxial approximation has been developed by Zhao and Gang (Zhao & Gang, 2009). The paper had some print errors, so here we show the correct formulation of the nonparaxial MGB. As below, we abbreviate "the nonparaxial MGB model" to "the MGB model".
A uniform normal pressure on the transducer surface can be expanded as the superposition of Gaussian beams (Wen & Breazeale, 1988). For a rectangular element with length 2a 1 and width 2a 2 in the x 1 -a n dx 2 -directions, respectively, the pressure over the transducer surface can be expressed as P 0 = P 0 |y 1 |≤a 1 , |y 2 |≤a 2 , y 3 = 0 0o t h e r w i s e = P 0 10 ∑ k 10 ∑ l A k A l exp(−B k y 2 1 /a 2 1 − B l y 2 2 /a 2 2 ) where A and B are ten complex coefficients obtained with an optimization method. Now the coordinate of wave field to be calculated is defined as . In order to obtain more accurate solution beyond the paraxial approximation region, the distance factor r is approximated as is not sensitive to a small element, so it can be substituted by K p (θ o ) a n dw em o v ei tt oo u t s i d eo ft h ei n t e g r a l . H e r e θ o is the angle between the vector (x − o) and the surface normal of transducer. Introducing 68 Ultrasonic Waves www.intechopen.com

3-D Modelings of an Ultrasonic Phased Array Transducer and Its Radiation Properties in Solid 11
Equations (23) and (24) into (14), we can obtain: Using the known integral formula the MGB model for a rectangular element can be obtained as where D 1 = k p a 2 1 /2 andD 2 = k p a 2 2 /2. To predict the wave field in solid from a rectangular element, we only perform the summation of 10 × 10 in Equation (27). Therefore fast calculation is possible but it is noted that the MGB is not applicable to complicated element shapes.
The total beam field of a phased array transducer can be calculated by simple superposition of individual wave fields with the corresponding time delay ∆τ nm . The P-wave displacement field can be calculated by where x is the position vector of the global coordinate. In Equation (27), the origin of x ′ is always set on the center of each element. The total wave field at the global coordinate can be obtained by adding up the results of Equation (27) in consideration of the delay and location of each element.

Numerical comparison of the RSNI and MGB model
Here the wave fields from an array transducer in stainless steel (c p =5800m/s, c s =3200m/s, ρ = 7800kg/m 3 ) are calculated with both the RSNI and the MGB model. The array transducer used in numerical comparisons is a linear phased array with 24 elements and center frequency f =3.0MHz. The element pitch p is 1.0mm and the gap g is 0.1mm. We show the magnitude of displacement |u| which is normalized by multiplying ρc 2 p /P 0 . In order to investigate the applicability of the MGB model to the UT modeling, the magnitudes of displacement fields are compared in both methods. Here we introduce the Rayleigh distance * : Fig.10(a) shows the on-axis (x 1 =0, x 2 =0) displacement when we vary F 3 =R, R/4 and R/8 with keeping the steering angle 0 • . From Fig.10(a), the result calculated by the MGB model can conform to the RSNI results. Similar behaviors are given in the off-axis results as shown in Fig.10(b). Fig.10(b) shows displacements on the x 1 -direction on x 3 =R, R/2 and R/4 in the case of F 3 =R, R/2 and R/4, respectively. It is understood that the MGB model is an effective tool in simulating the beam field radiated from phased arrays over a wide range of steering angle.

Visualization of 3-D wave field from a matrix array transducer
The following numerical examples are carried out with the MGB model because the MGB model can keep good accuracies in both the on-axis and off-axis regions. Here we show a visualization of 3-D wave field from a matrix array transducer. It is assumed that the array transducer has parameters as p x =p y =1.0mm and N × M=16 × 16=256. The 3-D wave fields in the cases of f =2.5 and 5.0MHz when we set the focal point to (F 1 ,F 2 ,F 3 )=(10mm,10mm,20mm) are illustrated in Fig.11. The main lobe only appears in the case of f =2.5MHz, however not only a main lobe but also grating lobes are generated in the case of f =5.0MHz. According to Equations (10) and (11), the directions of the grating lobe in f =5.0MHz are estimated to (θ ′ , φ ′ )=(58.8 • , −61.5 • ) and (58.8 • , 151.5 • ), on the other hand the grating lobe in f =2.5MHz is nonexistent. These estimation results show good agreements with the visualization results in Fig.11. * Note that it is well known that the actual focal length is shorter than the desired (designed) one and the focal point is never beyond the Rayleigh distance (Schmerr, 1998

Properties of linear and matrix array transducers
Here we investigate properties of both linear and matrix array transducers for the beam focusing and the beam steering.

Beam focusing
First, the property of beam focusing in the case of a linear array transducer is investigated. As a linear array transducer, the specification of p=1.0mm, g=0.1mm, N=24 and W/A=0.3 is assumed. The Rayleigh distance R of this transducer is 98.48mm in the case of f =4.0MHz. In Fig.12(a), the actual focal lengths F A are plotted in varying the desired focal length F 3 .T h e vertical axis in Fig.12(a) shows the parameter α(= F A /F 3 ) which represents the efficiency of the beam focus. When the parameter of α comes close to 1, it shows that we can transmit ultrasonic beam to the intended position. The parameter γ(= R/F 3 ) represents the closeness of focal length for the Rayleigh distance. From Fig.12(a), the focal length F A in the range of small γ is much shorter than the desired length F 3 .T h ep e a ko fα appears at γ = 4, and the distance between F A and F 3 increases when γ becomes large more than γ = 4.
Next let us look at the property of beam focusing for a matrix array transducer. Here a matrix array transducer with p x = p y =1.0mm, g x = g y =0.1mm, N=M=24 and A y /A x =1.0 is modeled. From Fig.12(b), the focal length F A comes closer and closer to F 3 when γ becomes larger and larger. Through the all range of γ, the difference between F A and F 3 is smaller than the case of the linear array transducer. Also ultrasonic beam can be focused efficiently even at the region near the surface of array transducer.

Beam steering
Figure 13(a) shows the deviation between the desired and actual focal points in a linear array transducer. Here the Rayleigh distance is R=98.48mm and we keep constant focal length |F |=R/4=24.62mm through all steering angles. The actual focus locations are less than the desired focal length and the deviation of the actual focus from desired focal position also increases as the steering angle increases. And Fig.13(b) shows maximal amplitude of normalized displacement and actual focal lengths F A versus the steering angle. These simulation results point out that the peak-to-peak amplitude will decrease as the steering angle increases. These changes will reduce the angular sensitivity of the phased array transducer.   In contrast, the behavior of the matrix array transducer is different from the linear array transducer. Figure 14(a) shows that the deviation distance between the desired and actual focal points is almost constant as the steering angle increases. However the peak-to-peak amplitude becomes lower as the steering angle increases.
From these results, the effective steering angle for the linear array transducer is up to approx. 40 • . Although the beam steering for the matrix array transducer is possible in higher angle than the linear array transducer, we suggest that the effective steering angle is up to approx. 55 • from the standpoint of the 6dB down of the peak-to-peak amplitude.

Modeling of phased array transducer in time domain
We here introduce a time domain simulation tool to predict the wave propagation from the array transducer. The method is a combined technique of the finite integration technique (FIT) (Marklein, 1998) and an image-based modeling (Terada et al., 1997). In the FIT, the finite difference equations of the stress and the particle velocity are derived from integral forms of the governing equations. Computational grids of these quantities are arranged in a staggered configuration in space, and the finite difference equations can be solved by marching time steps in the leap-frog manner. The FIT has an advantage that it can treat different boundary conditions without difficulty. This is essential to model the ultrasonic wave propagation in  heterogeneous material (Schubert & B. Köhler, 2001). In order to perform simulations of the UT accurately, a realistic shape data of the target is required. Here the image-based modeling is adopted as a pre-processing of the FIT. Using the image-base modeling, we can make geometries of targets directly from a digital image such as X-ray photograph, captured curve data of surface, CAD data, etc. Here we describe the 3-D image-based FIT (Nakahata et al., 2011) and show a numerical example of the phased array UT.

3-D finite integration technique (FIT)
We show the formulation of the FIT and the calculation flow as below. We consider the Cartesian coordinates (x 1 , x 2 , x 3 ). The governing equations of elastic waves are the Cauchy equation of motion and the equation of deformation rate. These equations are given in integral forms for a finite volume V with the surface S by where v is the particle velocity vector, τ is the stress second rank tensor, ρ is the density, n is the outward normal vector on the surface S,andf is the body force vector. In Equation (30), C is the stiffness tensor of rank four. In the case of isotropic materials, C can be written as:   in terms of the two Lamé constants, λ and µ. In Equation (32), δ ij is the Kronecker delta tensor. P and S wave velocities are given as: In the case of an acoustic problem, we can use above equations but have to set τ ij = 0 (i = j) in Equation (31) and µ = 0 in Equation (32).
ii integration cell Fig. 16. τ ii -integration cell (left figure) and v 1 -integration cell (right figure). All material parameters are defined in the τ ii -integration cell.
We consider a discretization form of Equation (30) using a cube as an integral volume V in Fig.15. Assuming that τ 11 is constant in V,wehavė whereτ 11 = ∂τ 11 /∂t and ∆x is the length on a side of V. In Equation (34), superscript () expresses positions of physical quantities in the integral volume as shown in Fig.16. Similarly a discretization of Equation (31) becomes where we let f = 0. Since the density ρ is given in the τ ii -integral volume as shown in Fig.16, the average value ρ =(ρ (L) + ρ (R) )/2 is used in Equation (35).
In the time domain, stress components τ are allocated at half-time steps, while the velocities v are at full-time steps. The following time discretization yields an explicit leap-frog scheme: where ∆t i st h et i m ei n t e r v a la n dt h es u p e r s c r i p tz denotes integer number of the time step. Therefore the FIT repeats the operations of Equations (36) and (37) by means of adequate initial and boundary conditions. A specific stability condition and adequate spatial resolution must be fulfilled to calculate the FIT accurately (Nakahata et al., 2011;Schubert & B. Köhler, 2001).

Image based modeling
We show a procedure of the 3-D image-based modeling with a 3-D curve measurement system based on the coded pattern projection method (Nakahata et al., 2012). The system consists of a projector that flashes narrow bands of black and white light onto the model's face and a camera that captures the pattern of light. The distortion pattern of light leads to the reconstruction of the surface of model. As shown in Fig.17, a 3-D shape of a welded T-joint is reconstructed from surface profiles which are measured at multiple angles. After unnecessary parts for calculation are trimmed away, the rest part is discretized into a voxel data. Here voxel is "volumetric pixel", representing a value on a regular grid in a 3-D space. The voxel data is directly used as the computational cell in the FIT, therefore the voxel size is set to be equivalent to the cell size of the FIT.

3-D numerical simulation
A 3-D simulation of ultrasonic wave propagation is demonstrated using the model in Fig.17.
Here an artificial defect which is assumed to be a lack of weld penetration is introduced at the inside of the welding part(see in Fig.18). The diameter of the defect is approx. 3mm and thickness is 0.1mm. Material parameters of the welded T-joint are set as c P = 5800m/s, c S = 3100m/s and ρ = 7800kg/m 3 . The linear array transducer (total element number N=24, the pitch p=0.6m, the element width W=7mm and the gap g=0.1mm) is located on the top of the model. In the simulation, we choose ∆x = 0.02mm, ∆t=1.0 nano-second (ns) and the total time step is 6500. Total number of the voxel is about 1 billion (1000 millions). The P-wave with the center frequency 4.0MHz is transmitted from the array transducer toward the defect. Here the ultrasonic wave is generated into solid by giving a time-dependent normal stress at the surface of the transducer. Figure 19 shows the isosurface of absolute value of the displacement vector |u(= vdt)| and Fig.20 shows |u| at a cross section of the model. These results show that ultrasonic wave propagates toward the artificial defect and scattered waves are generated. After some time steps, the echo from the defect arrives at the array transducer. In the general array UT, the echoes at the receiving process are also synthesized with the appropriate delay time in the same way to the transmitting process. Figure 21 shows the synthesized echo at receiving process, and waveform around 5 µs is the wave component from the defect.

Summary
For the reliable application of phased array techniques to NDT, it is essential to have thorough understanding on the properties of the phased array transducer and characteristics of radiation beam in solid. Here we show principle of electronic scanning with phased array transducers and an effective setting of the delay time to avoid grating and side lobes. In order to predict the wave field due to the array transducer, 3-D calculation tools in both frequency and time domains are introduced. From the simulation results with the MGB model in the frequency domain, characteristics about beam steering and focusing are compared between the linear and matrix array transducers. As the time domain calculation tool, the 3-D image-based FIT modeling of the UT for a welded T-joint is demonstrated.