# Family of Quantum Sources for Improving Near Field Accuracy in Transducer Modeling by the Distributed Point Source Method

Dominique Placko, Thierry Bore, Tribikram Kundu

► To cite this version:

Dominique Placko, Thierry Bore, Tribikram Kundu. Family of Quantum Sources for Improving Near Field Accuracy in Transducer Modeling by the Distributed Point Source Method. Applied Sciences, 2016, 6 (10), 10.3390/app6100302. [https://doi.org/10.3390/app6100302](https://doi.org/10.3390/app6100302)

HAL Id: hal-01793324

https://hal.science/hal-01793324

Submitted on 29 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Family of Quantum Sources for Improving Near Field Accuracy in Transducer Modeling by the Distributed Point Source Method

Dominique Placko<fnref n="1" />, Thierry Bore<fnref n="1" />,<fnref n="2" />,* and Tribikram Kundu<fnref n="3" />

Academic Editor: César M. A. Vasques

Received: 14 July 2016; Accepted: 10 October 2016; Published: 18 October 2016

**Abstract:** The distributed point source method, or DPSM, developed in the last decade has been used for solving various engineering problems—such as elastic and electromagnetic wave propagation, electrostatic, and fluid flow problems. Based on a semi-analytical formulation, the DPSM solution is generally built by superimposing the point source solutions or Green’s functions. However, the DPSM solution can be also obtained by superimposing elemental solutions of volume sources having some source density called the equivalent source density (ESD). In earlier works mostly point sources were used. In this paper the DPSM formulation is modified to introduce a new kind of ESD, replacing the classical single point source by a family of point sources that are referred to as quantum sources. The proposed formulation with these quantum sources do not change the dimension of the global matrix to be inverted to solve the problem when compared with the classical point source-based DPSM formulation. To assess the performance of this new formulation, the ultrasonic field generated by a circular planar transducer was compared with the classical DPSM formulation and analytical solution. The results show a significant improvement in the near field computation.

**Keywords:** acoustic source modeling; numerical method; distributed point source method

## 1. Introduction

Modeling of ultrasonic and sonic fields generated by a planar transducer of finite dimension is a fundamental problem whose solution is available in textbooks in this area [1–4]. A good review of the earlier developments of the ultrasonic field modeling in front of a planar transducer can be found in [5]. A list of more recent developments is given by Sha et al. [6]. The pressure field in front of a planar transducer in a homogeneous isotropic fluid has been computed both in the time domain [5,7,8] and in the frequency domain [9–14]. In the frequency domain analysis, which is more popular, the transducers are assumed to be vibrating with constant amplitudes at certain frequencies, and the pressure fields in front of the transducers are computed. Most of these models give the steady state response of the transducers. Recently phased array transducers have been modeled as well [15–20].

In all modeling approaches that are based on the Huygens-Fresnel superposition principle and Rayleigh-Sommerfeld Integral (RSI) the strengths of the point sources modeling the transducer are assumed to be known. For type-I RSI the pressure distribution, and for type II RSI, the velocity distribution, on the transducer surfaces are assumed to be known. The distributed point source method (DPSM) [21–23] has been developed relatively recently for solving electrostatic, electromagnetic,

<footnote n="1">¹ Ecole Normale Superieure, SATIE, 61, Avenue Du President Wilson, CEDEX F-94235 Cachan, France; dominique.placko@satie.ens-cachan.fr</footnote>

<footnote n="2">² School of Civil Engineering, Advanced Engineering Building, University of Queensland, St Lucia 4072, QLD, Australia</footnote>

<footnote n="3">³ Department of Civil Engineering and Engineering Mechanics, University of Arizona, 1209 E 2nd Street, Tucson, AZ 85721, USA; tkundu@email.arizona.edu</footnote>

<footnote>* Correspondence: bore.thierry@gmail.com; Tel.: +61-7-3365-8395</footnote>
and ultrasonic wave propagation problems. DPSM was first proposed by Placko and his co-workers for solving electrostatic and magnetostatic problems [21]. Later it was developed further by Placko and Kundu [22] and Kundu et al. [23] for solving ultrasonic problems. In ultrasonic modeling the source strengths can be treated as unknowns and obtained from the transducer surface boundary conditions as velocity or pressure. Alternately, all point sources can be assumed to have the same strength in order to model the uniform velocity condition at the transducer surface. In recent years DPSM has been extended to solve high frequency fluid flow problems [24] and fluid mechanics problems involving flow around an obstacle [25].

DPSM showed a significant savings in computational time compared to the Finite Elements Method (FEM) FEM technique for modeling ultrasonic problems [23,24,26] or electromagnetic problems [27]. A three-dimensional (3D) wave field modeling problem in a homogeneous fluid in front of a square transducer was solved in 2 min by DPSM while COMSOL Multi-Physics FEM code took 35 h to solve the identical problem on the same computer [23]. Jarvis and Cegla [26] used DPSM and FEM for solving ultrasonic wave reflections from a rough surface and concluded, "[FEM solution is] two orders of magnitude slower than the equivalent DPSM simulation on the same machine".

In the above works on DPSM, point sources were uniformly distributed near the boundaries, interfaces, and front faces of the transducer (or emitter). The source strengths were then calculated by solving a system of linear equations that are obtained by satisfying boundary and continuity conditions for the problem geometry. In the near field, very close to the point sources, the computed field satisfies the boundary conditions only at the discrete points where the boundary conditions are specified. In between those points the boundary conditions are generally not satisfied. To overcome this shortcoming the use of a scaling factor was proposed [23]. In [23] Kundu et al. showed that the computed values in the far field (away from the transducer face) could be accurately modeled after incorporating the scaling factor. After considering this scaling factor for both near and far field computations, it was observed that at sections very close to the transducer face the average value of the computed field matched well with the true solution, although for individual points the deviation from the true solution could still be significant. To overcome this shortcoming Rahani and Kundu [28] tried DPSM with equivalent sources densities instead of the unique point source. They proposed two modifications—(1) G-DPSM or Gaussian DPSM replacing the single point source by a family of sources whose strength variations follow a Gaussian profile, and (2) ESM (Element Source Method): replacing discrete point sources by continuously varying element sources. Although these two modifications avoided the near field inaccuracy problem both approaches had their limitations. Results generated by G-DPSM were strongly affected by the Gaussian distribution function parameters and the ESM formulation was relatively complex. Thus, one major advantage of the DPSM formulation, its simplicity, was compromised by ESM. The concept of equivalent source density, including elemental point sources, linear, surface, or volume repartition, was presented and advocated by Placko et al. [29].

The objective of this paper is to propose a new ESD (equivalent source density) approach which is much simpler and can improve the near field computation accuracy without increasing the matrix size that is to be inverted for obtaining the solution. In this solution, the single point source is replaced by a source family. Point sources in a family are denoted as quantum sources to distinguish them from the regular point sources used in the classical DPSM formulation. A transducer can be modeled by *M* number of regular point sources in classical DPSM formulation or *M* × *F* quantum sources in the modified formulation. If every one of the *M* point sources are replaced by *F* number of quantum sources then the total number of quantum sources is equal to *M* × *F*. The quantum source strengths are set such that a given radiation pattern from the quantum source family is obtained.

G-DPSM, proposed by Rahani and Kundu [28], assumed a Gaussian distribution of source strength variations in a source family and their final results were strongly affected by the Gaussian distribution function parameters. However, in the technique proposed here, the strength variation of the quantum sources in a family is assumed to be unknown, but the radiation pattern generated by a source family is considered to follow a specific pattern—Gaussian distribution is one possibility,
but not necessarily the only possibility. However, for not assuming a Gaussian distribution of the source strengths in a family significantly increases the number of unknown source strengths. It is illustrated in the Formulations section how this improvement can be achieved without increasing the final matrix size that is to be inverted. Thus, this modification can improve the computational accuracy without significantly increasing the computational time.

## 2. Formulations

The classical and the proposed DPSM formulations using family ESD are presented in this section. Although the classical DPSM formulation using unique point sources is available in the literature—in various publications mentioned above—for the benefit of the readers the basic steps of classical DPSM formulation and its current limitations are briefly reviewed first in the following subsection, before introducing the new formulation.

### 2.1. Classical DPSM Formulation: ESD is a Single Point Source

To generate the acoustic field in the fluid medium in front of an ultrasonic transducer or the electric field in front of an emitter by the DPSM technique, a number of point sources are placed slightly behind the transducer (or emitter) face as shown in Figure 1. This figure shows M number of source points placed at the centers of M spheres of radius $r_s$. Spheres touch the front face of the transducer. Thus, the point sources are located at a distance of $r_s$ behind the transducer face to avoid the singularity. The DPSM formulation for ultrasonic wave field modeling is presented below. However, for electric field modeling the derivations are identical; the parameter names are simply changed from pressure, velocity, etc., to electric field, magnetic field, etc.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** M point sources are placed slightly behind the transducer face and the ultrasonic field is computed at N target points in front of the transducer.</figcaption></figure>

From [3] the pressure and velocity fields in a fluid at point x at a distance $r_m$ from the m-th point source of strength Am can be written as:

The pressure field:

$$p_m(\mathbf{x}) = p_m(r_m) = A_m \frac{\exp(ik_f r_m)}{r_m} \qquad (1)$$

The velocity in the radial direction, at a distance r from the m-th point source:

$$v_m(r) = \frac{A_m \exp(ik_f r)}{i\omega\rho r} \left( ik_f - \frac{1}{r} \right) \qquad (2)$$
Therefore, the z-component of the velocity is:

$$v_{zm}(r) = \frac{A_m}{i\omega\rho} \frac{\exp(ik_f r)}{r^2} \left(ik_f - \frac{1}{r}\right) \quad (3)$$

If there are *M* point sources distributed over the transducer surface, as shown in Figure 1, then the total pressure and the z-direction velocity values at point *x* are given by:

$$p(x) = \sum_{m=1}^{M} p_m(r_m) = \sum_{m=1}^{M} A_m \frac{\exp(ik_f r_m)}{r_m} \quad (4)$$

$$v_z = \sum_{m=1}^{M} v_{zm}(r_m) = \sum_{m=1}^{M} \frac{A_m}{i\omega\rho} \frac{z_m \exp(ik_f r_m)}{r_m^2} \left(ik_f - \frac{1}{r_m}\right) \quad (5)$$

where $z_m$ is the distance of point *x* in the z-direction (perpendicular to the transducer face) measured from the *m*-th point source.

If the transducer surface velocity in the z direction is given by $v_0$ then for all *x* values on the transducer surface the velocity should be equal to $v_0$. Therefore:

$$v_z(x_{\text{transducer}}) = \sum_{m=1}^{M} \frac{A_m}{i\omega\rho} \frac{z_m \exp(ik_f r_m)}{r_m^2} \left(ik_f - \frac{1}{r_m}\right) = v_0 \quad (6)$$

By satisfying Equation (6) at *M* points on the transducer surface, where the spheres touch the transducer face it is possible to obtain a system of *M* linear equations with *M* unknowns ($A_1, A_2, A_3, \dots, A_m$). If the point sources are placed slightly behind the transducer face as shown in Figure 1 then $r_m$ of Equation (6) can never be equal to zero; thus, singular points are avoided.

If point *x* in Figure 1 coincides with the *n*-th target point and the velocity at the *n*-th target point or the observation point is to be computed then from Equation (6) one obtains:

$$v_z(x_n) = \sum_{m=1}^{M} \frac{A_m}{i\omega\rho} \frac{z_{nm} \exp(ik_f r_{nm})}{r_{nm}^2} \left(ik_f - \frac{1}{r_{nm}}\right) \quad (7)$$

In Equation (7) the first subscript *n* (of *z* and *r*) corresponds to the *n*-th target point and the second subscript *m* represents the *m*-th source point; therefore, $r_{nm}$ is the distance of the *n*-th target point from the *m*-th source point while $z_{nm}$ is the distance along the z-direction between the *n*-th target point and *m*-th source point. If the number of target points is made equal to the number of source points (both equal to *M*) and the target points are placed on the transducer surface then Equation (7) can be written in matrix form as:

$$V_S = M_{SS}A_S \quad (8)$$

where $V_S$ is the ($M \times 1$) vector of the velocity components at *M* number of target points on the transducer surface, and $A_S$ is the ($M \times 1$) vector containing the strengths of *M* number of point sources. $M_{SS}$ is the ($M \times M$) square matrix relating the two vectors $V_S$ and $A_S$:

$$\{V_S\}^T = \begin{bmatrix} v_{z1} & v_{z2} & \dots & v_{zM} \end{bmatrix} \quad (9)$$

where superscript *T* indicates the transpose and $v_{zm}$ represents the z-direction velocity at the *m*-th target point. If the transducer surface vibrates with a constant velocity amplitude $v_0$ then all elements of Equation (9) should have the same value $v_0$.
Vector $A_S$ of the source strengths is given by:

$$ \{A_S\}^T = \begin{bmatrix} A_1 & A_2 & A_3 & \dots & A_{M-2} & A_{M-1} & A_M \end{bmatrix} \qquad (10) $$

The matrix $M_{SS}$ is obtained from Equation (7):

$$ M_{SS} = \begin{bmatrix} m_{11} & m_{12} & \dots & m_{1M} \\ m_{21} & m_{22} & \dots & m_{2M} \\ \vdots & \vdots & \ddots & \vdots \\ m_{M1} & m_{M2} & \dots & m_{MM} \end{bmatrix} = \begin{bmatrix} f(z_{11}, r_{11}) & f(z_{12}, r_{12}) & \dots & f(z_{1M}, r_{1M}) \\ f(z_{21}, r_{21}) & f(z_{22}, r_{22}) & \dots & f(z_{2M}, r_{2M}) \\ \vdots & \vdots & \ddots & \vdots \\ f(z_{M1}, r_{M1}) & f(z_{M2}, r_{M2}) & \dots & f(z_{MM}, r_{MM}) \end{bmatrix}_{M \times M} \qquad (11) $$

where:

$$ m_{nm} = f(z_{nm}, r_{nm}) = \frac{z_{nm}\exp(ik_f r_{nm})}{i\omega\rho r_{nm}^2} \left( ik_f - \frac{1}{r_{nm}} \right) \qquad (12) $$

Equation (8) can be inverted to obtain the point source strengths:

$$ A_S = [M_{SS}]^{-1} V_S \qquad (13) $$

Instead of velocity, if pressure is specified on the transducer surface, then the source strength is obtained from the following equation:

$$ A_S = [Q_{SS}]^{-1} P_S \qquad (14) $$

where $P_S$ is the specified pressure field at the transducer surface and $Q$ matrix relates the pressure field and the source strength (obtained from Equation (4)).

It should be noted here that Equation (13) can also be rewritten as $V_S^* = M_{SS}^* A_S^*$, where $V_S^*$ is a $(3M \times 1)$ vector containing all three components of velocity at $M$ test points. For inversion of this equation, when needed, we may have used the concept of ESD, which are triplet sources used at every center point so that $A_S^*$ becomes a $(3M \times 1)$ vector containing all three strengths of the triplets at $M$ source center points. See [3] for more detailed derivation.

After evaluating the source strength vector $A_S$ from Equation (13) or (14) the pressure $p(x)$ or velocity vector $V(x)$ at any set of target points $x$ (on the transducer surface or away) can be obtained from the following equations:

$$ \begin{aligned} P_T &= Q_{TS}A_S \\ V_T &= M_{TS}A_S \\ V_T^* &= M_{TS}^*A_S \end{aligned} \qquad (15) $$

For a set of $N$ target points the vector $P_T$ contains $N$ pressure values at $N$ target points, so its size is $(N \times 1)$. $V_T$ is a $(N \times 1)$ vector containing z-components of velocity at $N$ target points, and $V_T^*$ is a $(3N \times 1)$ vector containing all three components of velocity at $N$ target points.

Therefore, the velocity field on the transducer surface can be computed at $N$ number of points where $N$ is much greater than $M$. It should be noted that only at $M$ points the boundary conditions are specified. The velocity variation on the transducer face shown in Figure 2 is obtained with $M = 464$ (the number of source points considered for modeling the transducer) and $N = 2500$ (the number of target points where velocity is computed on the plane of the transducer face). Velocity values are computed on a square surface that coincides with the transducer face. Some points on the square surface that are near its four corners are outside the transducer face and give zero displacement, as expected. However, on the transducer face the computed velocity is not uniform. It reaches a peak value of $v_0$ (in our case 1) at 464 points where small spheres of Figure 1 touch the transducer face but, in between these points, they are significantly smaller. Thus, the average velocity of the transducer face becomes less than $v_0$ giving rise to smaller pressure values from the DPSM analysis. When the DPSM results are multiplied by an appropriate scaling factor (obtained by dividing the given transducer
face velocity by the average velocity of the transducer surface) then the discrepancy in the far field disappears [3]. However, in the near field (on the transducer face and very close to the transducer face) the scaling factor multiplication cannot remove the oscillation patterns shown in Figure 2 although their average value matches with the prescribed velocity value of the transducer face.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Normal velocity on the transducer face computed by the DPSM (distributed point source method) technique when point sources, as shown in Figure 1, are used to model the transducer face.</figcaption></figure>

## 2.2. Formulation with ESD Families

In the proposed formulation every point source of Figure 1 is replaced by a family of point sources as shown in Figure 3. These point sources in the family are called quantum sources to distinguish them from the unique point sources in the classical DPSM formulation. If there are *F* number of quantum sources in every family then, from *M* number of original point sources, a total of (*M* × *F*) quantum sources are generated. If classical DPSM formulation as described above is applied to these (*M* × *F*) quantum sources then the matrix size becomes (*M* × *F*) × (*M* × *F*) (see Equations (13) and (14)).

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Quantum source clusters or source families replace the individual point sources in the new formulation for modeling the transducer.</figcaption></figure>
To reduce the matrix size and the number of unknown source strengths from (M × F) to M it is assumed that the relative source strength distribution within a family is the same for all M families. Then, F sources of the i-th family (i = 1, 2, ..., M) should have the following source strength distribution:

$$A_i = \lambda_i q = \lambda_i \begin{bmatrix} q_1 & q_2 & q_3 & \cdots & q_{(F-1)} & q_F \end{bmatrix}^T \quad (16)$$

The i-th family strength is denoted by $\lambda_i$ and its distribution is denoted by the vector q. The distribution vector q is the same for all M families but the strength $\lambda_i$ can be different for different families. Then the velocity field at M target points generated by (M × F) quantum sources can be written as:

$$\left\{ \begin{matrix} v_{z1} \\ v_{z2} \\ v_{z3} \\ \vdots \\ v_{zM} \end{matrix} \right\} = \left[ \begin{matrix} m_{11}^1 & \cdots & m_{1F}^1 \\ m_{21}^1 & \cdots & m_{2F}^1 \\ m_{31}^1 & \cdots & m_{3F}^1 \\ \vdots & & \vdots \\ m_{M1}^1 & \cdots & m_{MF}^1 \end{matrix} \begin{matrix} m_{11}^2 & \cdots & m_{1F}^2 \\ m_{21}^2 & \cdots & m_{2F}^2 \\ m_{31}^2 & \cdots & m_{3F}^2 \\ \vdots & & \vdots \\ m_{M1}^2 & \cdots & m_{MF}^2 \end{matrix} \begin{matrix} m_{11}^M & \cdots & m_{1F}^M \\ m_{21}^M & \cdots & m_{2F}^M \\ m_{31}^M & \cdots & m_{3F}^M \\ \vdots & & \vdots \\ m_{M1}^M & \cdots & m_{MF}^M \end{matrix} \begin{matrix} q_1 \\ \vdots \\ q_F \\ \vdots \\ q_1 \\ \vdots \\ q_F \end{matrix} \right] \quad (17)$$

where $m_{jk}^i$ represents the velocity at the i-th target point generated by the k-th quantum source of unit amplitude at the j-th family. The above equation can be re-written as:

$$\left\{ \begin{matrix} v_{z1} \\ v_{z2} \\ v_{z3} \\ \vdots \\ v_{zM} \end{matrix} \right\} = \left[ \begin{matrix} (m_{11}^1 \quad \cdots \quad m_{1F}^1) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & (m_{11}^2 \quad \cdots \quad m_{1F}^2) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & \cdots & (m_{11}^M \quad \cdots \quad m_{1F}^M) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} \\ (m_{21}^1 \quad \cdots \quad m_{2F}^1) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & (m_{21}^2 \quad \cdots \quad m_{2F}^2) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & \cdots & (m_{21}^M \quad \cdots \quad m_{2F}^M) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} \\ (m_{31}^1 \quad \cdots \quad m_{3F}^1) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & (m_{31}^2 \quad \cdots \quad m_{3F}^2) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & \cdots & (m_{31}^M \quad \cdots \quad m_{3F}^M) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} \\ \vdots & \vdots & \ddots & \vdots \\ (m_{M1}^1 \quad \cdots \quad m_{MF}^1) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & (m_{M1}^2 \quad \cdots \quad m_{MF}^2) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} & \cdots & (m_{M1}^M \quad \cdots \quad m_{MF}^M) \begin{pmatrix} q_1 \\ \vdots \\ q_F \end{pmatrix} \end{matrix} \right] \left\{ \begin{matrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_M \end{matrix} \right\} \quad (18)$$

Or:

$$\left\{ \begin{matrix} v_{z1} \\ v_{z2} \\ v_{z3} \\ \vdots \\ v_{zM} \end{matrix} \right\} = \left[ \begin{matrix} M_{11} & M_{12} & \cdots & M_{1M} \\ M_{21} & M_{22} & \cdots & M_{2M} \\ M_{31} & M_{32} & \cdots & M_{3M} \\ \vdots & \vdots & \ddots & \vdots \\ M_{M1} & M_{M2} & \cdots & M_{MM} \end{matrix} \right] \left\{ \begin{matrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_M \end{matrix} \right\} \quad (19)$$

Therefore:

$$\left\{ \begin{matrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \vdots \\ \lambda_M \end{matrix} \right\} = \left[ \begin{matrix} M_{11} & M_{12} & \cdots & M_{1M} \\ M_{21} & M_{22} & \cdots & M_{2M} \\ M_{31} & M_{32} & \cdots & M_{3M} \\ \vdots & \vdots & \ddots & \vdots \\ M_{M1} & M_{M2} & \cdots & M_{MM} \end{matrix} \right]^{-1} \left\{ \begin{matrix} v_{z1} \\ v_{z2} \\ v_{z3} \\ \vdots \\ v_{zM} \end{matrix} \right\} \quad (20)$$
It should be noted that both Equations (13) and (20) require inversion of the $M \times M$ matrix. However, for classical DPSM formulation (Equation (13)) the number of point sources is *M* while for the new formulation (Equation (20)) the number of quantum sources is (*M* × *F*).

To obtain the relative strength variation vector **q** = [ q₁ q₂ q₃ ... qᵢⱼ ... qᵣⱼⱼ ]ᵀ (see Equation (16)) of the quantum sources in a family the radiation pattern generated by the quantum sources will have to be assumed. The radiation pattern describes the velocity or pressure distribution generated by a family of quantum sources for ultrasonic field modeling, or electric field distribution for electrostatic problems.

In the classical DPSM formulation Green's function for a point source determines the radiation pattern. Typically the wavefront generated by a point source is a sphere in an isotropic medium. However, for a family of quantum sources any radiation pattern can be achieved by properly adjusting the source strength distribution vector **q** in the family. We may notice that, for the far field, all of the terms of the line matrix ( m₁₁ ... m₁ⱼ ... mᵢⱼ ... mᵣⱼⱼ ) of Equation (18) tends to the same value, which corresponds to the *m*₁₁ term of Equation (12). Then, the family is seen as a cluster of sources equivalent to a unique point source of strength equal to the sum of [ q₁ q₂ q₃ ... qᵢⱼ ... qᵣⱼⱼ ].

### 3. Results

Numerical results are generated for a Gaussian radiation pattern, as shown in Figure 4, generated by every family. Quantum source strengths within every family are obtained for this radiation pattern following the classical DPSM approach and ignoring the presence of other families. Then the pressure field in the fluid in front of the uniformly vibrating circular transducer with a radius of 6 mm is computed by considering all families, as shown in Figure 3. For clarity, 50 families are shown in Figure 3. However, during actual computation 464 families were considered. Every family is composed of 30 quantum sources. The transducer vibration frequency is 1 MHz. Fluid properties considered for this simulation are those of water—density (ρ) = 1000 kg/m³ and P-wave speed (c<sub>f</sub>) = 1480 m/s.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Normal velocity variation (or radiation pattern) from one source family.</figcaption></figure>

The normal velocity component is first computed at 2500 points on a square surface that coincides with the transducer surface. The computed results with the new formulation is shown in Figure 5a and with the old formulation (classical DPSM) is shown in Figure 5b. Clearly, the oscillations are much less (almost disappearing) in Figure 5a. In both formulations the final matrix that was inverted had dimensions of 464 × 464. In the new formulation 464 clusters gave this matrix dimension, while in the old formulation 464 point sources generated this matrix. To improve the accuracy of the classical DPSM computation, the number of point sources and, hence, the matrix dimension must be increased while, in the new formulation, the number of quantum sources can be increased without increasing the number of clusters and, therefore, the matrix dimension.
<figure><img src="image_6.png" /><figcaption>Figure 5. The normal velocity on the transducer face computed by the DPSM technique when the transducer face is modeled by (a) source clusters, as shown in Figure 3 (new formulation); and (b) individual point sources as shown in Figure 1 (old formulation).</figcaption></figure>

Figure 6 shows the variation of the normal velocity along the diameter (y-axis going through the center of the transducer face). The continuous line shows the field computed by the new formulation while the dotted line is obtained from the classical DPSM formulation. Nine points where the two curves match correspond to the boundary points where the velocity is specified. These nine points out of 464 points appear on the diameter along which the velocity values are computed. Classical DPSM gives correct values only at those nine points and deviates significantly from the correct value of 1 at other points. The results generated by the new formulation are much closer to the true values and match with the true values at several other points where the boundary conditions are not specified. To get a quantitative measure of the difference between these curves the error is computed using the following equation:

$$ Error = 100 \times \sqrt{\frac{\sum_{i} \{V(i)_{\text{true}} - V(i)_{\text{computed}}\}^2}{\sum_{i} V(i)_{\text{true}}^2}} \quad (21) $$

where $V(i)_{\text{true}}$ is the true value of the velocity amplitude at the i-th point, which is 1 on the transducer surface and zero outside the transducer, as shown in Figure 6. The true value is denoted as the objective function (dashed-dotted line) in this figure. $V(i)_{\text{computed}}$ are the computed values (continuous and dotted lines of Figure 6). The computed error is 34.8% for classical DPSM and 3.05% for proposed DPSM with the family ESD technique.

<figure><img src="image_7.png" /><figcaption>Figure 6. The normal velocity variation on the transducer surface: dashed-dotted line—true value; solid line—computed results using source clusters, as shown in Figure 3; and dashed line—computed by superimposing individual point source contributions, as shown in Figure 1.</figcaption></figure>
From the symmetry of the problem it is obvious that the fluid velocity parallel to the transducer surface in the x-direction should be zero along the y-axis going through the transducer center. This velocity component on the transducer surface is plotted in Figure 7. It is close to zero for the new formulation. However, for the classical formulation, strong velocity values are observed near the periphery of the transducer. Both formulations give zero velocity near the central axis of the transducer.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** The X-component of the velocity variation on the transducer surface: solid line—computed results using source clusters, as shown in Figure 3; and dashed line—computed by superimposing individual point source contributions, as shown in Figure 1.</figcaption></figure>

From Figures 6 and 7 it is clear that the new formulation significantly improves the field computed in the fluid adjacent to the transducer face. To investigate how the computed fields away from the transducer face are affected by these two formulations the normal velocity components are plotted along the y-axis for four different distances (z = 0, r<sub>s</sub>/3, 2r<sub>s</sub>/3 and r<sub>s</sub>) from the transducer surface. These figures are shown in Figure 8. It should be noted that r<sub>s</sub> is the radius of the source sphere, as shown in Figure 1. Therefore, r<sub>s</sub> is also the distance of the point source from the transducer surface. Figure 8 clearly shows that as the observation points move away from the transducer face by a distance equal to r<sub>s</sub>, or when they are located at a distance 2r<sub>s</sub> from the point sources, then the oscillations in the computed field disappear for the classical DPSM also. However, the curves generated by the classical DPSM and the family ESD DPSM do not match, although their shapes appear to be similar. The classical DPSM results need to be multiplied by a scaling factor to match with the new results.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** The normal velocity variation on the transducer surface (bottom figure) and near it (the other three figures), 'z' is the distance from the transducer face; solid line—computed results using source clusters, as shown in Figure 3; and dashed line—computed by superimposing individual point source contributions, as shown in Figure 1.</figcaption></figure>
Classical and family ESD DPSM results are then compared with the analytical solution in Figure 9. For a circular transducer of radius of the analytical expression of the pressure field along the central axis is given by [23]:

$$p(z, \omega) = \rho v_0 c_f \left\{ \exp(ik_f z) - \exp(ik_f \sqrt{z^2 + a^2}) \right\} \quad (22)$$

where *p* is the fluid pressure at a distance *z* from the transducer face vibrating at frequency *ω* (rad/s) in a perfect fluid of density *ρ*. The normal velocity component (*v*<sub>z</sub>) at the transducer surface is *v*<sub>0</sub>. The wave number *k*<sub>*f*</sub> is equal to *ω*/*c*<sub>*f*</sub>, where *c*<sub>*f*</sub> is the acoustic wave speed in the fluid. The dashed-dotted line of Figure 9 shows the analytical solution, the continuous line is the family DPSM solution, and the dotted line is the classical DPSM solution. Clearly, the family DPSM solution is closer to the analytical solution (error is 5.57%) than the classical DPSM solution (error is 32.7%), but neither solution matches with the analytical solution. This is because we arbitrarily assumed the Gaussian radiation pattern for each family. This assumption might be better than the spherical radiation pattern generated by a point source in classical DPSM modeling but it is still not necessarily the actual radiation pattern needed for perfect matching. However, both curves can be matched with the analytical solution simply by multiplying them with appropriate scaling factors.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Pressure variation along the central axis in front of the transducer: dashed-dotted line—theoretical value; solid line—computed results using source clusters, as shown in Figure 3; and dashed line—computed by superimposing individual point source contributions, as shown in Figure 1.</figcaption></figure>

The sensitivity of the family DPSM-computed results on the assumed radiation pattern of the individual families is investigated in Figures 10 and 11. The standard deviation (σ) of the Gaussian radiation profile is changed from 1.3r<sub>s</sub> to 1.7r<sub>s</sub> with an interval of 0.2r<sub>s</sub> and the normal velocity profiles on the transducer surface along a diameter are plotted in Figure 10. Clearly, the results are sensitive to the assumed radiation pattern. However, all three curves generated by the family DPSM formulation are much closer to the true solution (v<sub>z</sub> = 1) than the classical DPSM-generated results. It should be mentioned here that Figures 5–9 were generated for σ = 1.5r<sub>s</sub>.

Finally, Figure 11 shows the pressure field variations along the central axis generated by classical DPSM and family DPSM with same three values of σ as considered in Figure 10. Note that for some radiation pattern the family DPSM-computed field matches very well with the analytical solution. However, all of these curves can be brought very close to the true solution by multiplying them with the appropriate scaling factor.
These results show that if we are interested in computing the field very close to the transducer surface ($z < r_s$) then the family DPSM formulation must be adopted because simply multiplying the computed field by a scaling factor cannot get rid of the oscillations in the computed field near the transducer. However, if such near field solutions are not needed, then the classical DPSM solution with a scaling factor generates accurate results.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** The normal velocity variation on the transducer surface: dashed line—true value; dotted line—computed by superimposing individual point source contributions, as shown in Figure 1; solid lines—computed results using source clusters, as shown in Figure 3, for three different radiation patterns.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Pressure variation along the central axis in front of the transducer: dashed line—theoretical value; dotted line—computed by superimposing individual point source contributions, as shown in Figure 1; solid lines—computed results using source clusters, as shown in Figure 3, for three different radiation patterns.</figcaption></figure>

4. Conclusions

A new formulation is presented to improve the accuracy of DPSM for the near field computation. The classical DPSM technique works very well for the far field computation but cannot accurately model the near field—the field on the transducer surface and very close to the transducer surface—because DPSM is based on the Huygens-Fresnel principle, assumes that the point of observation is not very close to the transducer face modeled by distributed point sources. The new formulation explores a new possibility of equivalent source density and makes the DPSM technique accurate in both the near field and the far field. This is achieved by replacing every point source in the classical
DPSM formulation by families of quantum sources. An individual family of sources can generate any radiation pattern, unlike an individual point source for which the radiation wave front can only be spherical in an isotropic medium. The new formulation is developed in such a manner that it does not increase the dimension of the square matrix that needs to be inverted for the source strength calculation. Thus, the main advantage of DPSM, that it requires low computational time compared to FEM, is not compromised. In both FEM and classical DPSM formulations the standard approach for improving accuracy is to increase the mesh size for FEM or the number of point sources for DPSM. Both of these approaches significantly increase the computational time and is avoided here. Numerical examples are presented demonstrating the improvement in the computational accuracy with the new technique.

**Acknowledgments:** The third author acknowledges the financial support from Ecole Normale Supérieure de Cachan, France towards his visit as an Invited Professor to conduct this collaborative research.

**Author Contributions:** All three authors contributed to the paper in terms of conceiving the research idea, generating the numerical results and writing the paper. However, Dominique Placko took the leading role in conceiving the idea, Thierry Bore was leading in generating the numerical results and Tribikram Kundu took the leading role in writing the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Rayleigh, L. *Theory of Sound*; Dover Publications: New York, NY, USA, 1965; Volume 2, pp. 162–169.

2. Morse, P.M.; Ingard, U.K. *Theoretical Acoustics*; McGraw-Hill: New York, NY, USA, 1968.

3. Schmerr, L.W. *Fundamental of Ultrasonic Nondestructive Evaluation—A Modeling Approach*; Plenum Press: New York, NY, USA, 1998.

4. Kundu, T. (Ed.) *Ultrasonic Nondestructive Evaluation: Engineering and Biological Material Characterization*; CRC Press: Boca Raton, FL, USA, 2004.

5. Harris, G.R. Review of transient field theory for a baffled planar piston. *J. Acoust. Soc. Am.* **1981**, *70*, 10–20.
[CrossRef]

6. Sha, K.; Yang, J.; Gan, W.S. A complex virtual source approach for calculating the diffraction beam field generated by a rectangular planar source. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **2003**, *50*, 890–895.
[PubMed]

7. Stepanishen, P.R. Transient radiation from piston in an infinite planar baffle. *J. Acoust. Soc. Am.* **1971**, *49*, 1627–1638.
[CrossRef]

8. Jensen, J.A.; Svendsen, N.B. Calculation of pressure fields from arbitrary shaped, apodized, and excited ultrasound transducers. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **1992**, *39*, 262–267.
[CrossRef]
[PubMed]

9. Ingenito, F.; Cook, B.D. Theoretical investigation of the integrated optical effort produced by sound field radiated from plane piston transducers. *J. Acoust. Soc. Am.* **1969**, *45*, 572–577.
[CrossRef]

10. Lockwood, J.C.; Willette, J.G. High-speed method for computing the exact solution for the pressure variations in the near field of a baffled piston. *J. Acoust. Soc. Am.* **1973**, *53*, 735–741.
[CrossRef]

11. Scarano, G.; Denisenko, N.; Matteucci, M.; Pappalardo, M. A new approach to the derivation of the impulse response of a rectangular piston. *J. Acoust. Soc. Am.* **1985**, *78*, 1109–1113.
[CrossRef]

12. Hah, Z.G.; Sung, K.M. Effect of spatial sampling in the calculation of ultrasonic fields generated by piston transducers. *J. Acoust. Soc. Am.* **1992**, *92*, 3403–3408.
[CrossRef]

13. Wu, P.; Kazys, R.; Stepinski, T. Analysis of the numerically implemented angular spectrum approach based on the evaluation of two-dimensional acoustic fields. Part I. Errors due to the discrete Fourier transform and discretization. *J. Acoust. Soc. Am.* **1995**, *99*, 1139–1148.
[CrossRef]

14. Lerch, T.P.; Schmerr, L.W.; Sedov, A. Ultrasonic beam models: An edge element approach. *J. Acoust. Soc. Am.* **1998**, *104*, 1256–1265.
[CrossRef]

15. Azar, L.; Shi, Y.; Wooh, S.C. Beam focusing behavior of linear phased arrays. *NDT E Int.* **2000**, *33*, 189–198.
[CrossRef]

16. Buiochi, F.; Martinez, O.; Ullate, L.G.; Espinosa, F.M. 3D computational method to study the focal laws of transducer arrays for NDE applications. *Ultrasonics* **2004**, *42*, 871–876.
[CrossRef]
[PubMed]
17. Song, S.J.; Kim, C.H. Simulation of 3-D radiation beam patterns propagated through a planar interface from ultrasonic phased array transducers. *Ultrasonics* **2002**, *40*, 519–524. [CrossRef]

18. Ahmad, R.; Kundu, T.; Placko, D. Modeling of phased array transducers. *J. Acoust. Soc. Am.* **2005**, *117*, 1762–1776. [CrossRef] [PubMed]

19. Park, J.S.; Song, S.J.; Kim, H.J. Calculation of radiation beam field from phased array ultrasonic transducers using expanded Multi-Gaussian beam model. In *Advances in Safety and Structural Integrity 2005*; Book Series: Solid State Phenomena; Trans Tech Publications: Zürich, Switzerland, 2006; Volume 110, pp. 163–168.

20. Zhao, X.; Gang, T. A study on non paraxial beam model and ultrasonic measurement model of phased arrays. *Ultrasonics* **2008**, *49*, 126–130. [CrossRef] [PubMed]

21. Placko, D.; Liebeaux, N. Modèle à sources Réparties: Un concept pour le calcul du champ émis par les circuits magnétiques ouverts. In *Acte du Deuxième Colloque Interdisciplinaire en Instrumentation (C2I)*; Hermes: Paris, France, 2001; pp. 523–530.

22. Placko, D.; Kundu, T. Modeling of ultrasonic field by distributed point source method. In *Ultrasonic Nondestructive Evaluation: Engineering and Biological Material Characterization*; Kundu, T., Ed.; CRC Press: Boca Raton, FL, USA, 2004; pp. 143–202.

23. Kundu, T.; Placko, D.; Rahani, E.K.; Yanagita, T.; Dao, C.M. Ultrasonic field modeling: A comparison between analytical, semi-analytical and numerical techniques. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **2010**, *57*, 2795–2807. [CrossRef] [PubMed]

24. Wada, Y.; Kundu, T.; Nakamura, K. Mesh-free distributed point source method for modeling viscous fluid layer motion between disks vibrating at ultrasonic frequency. *J. Acoust. Soc. Am.* **2014**, *136*, 466–474. [CrossRef] [PubMed]

25. Placko, D.; Rivollet, A.; Bore, T. Procédé et Dispositif de Modélisation D'une Grandeur Physique d'un Écoulement Fluide (Process and set-up for Modeling a Physical Value in a Fluid Flow). French Patent Appl. 1454675, 23 May 2014.

26. Jarvis, A.J.C.; Cegla, F.B. Scattering of near normal incidence SH waves by sinusoidal and rough surfaces in 3-D: Comparison to the scalar wave approximation. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **2014**, *61*, 1179–1190. [CrossRef] [PubMed]

27. Bore, T.; Placko, D.; Taillade, F.; Sabouroux, P. Electromagnetic characterization of grouting materials of bridge post tensioned ducts for NDT using capacitive probe. *NDT E Int.* **2013**, *60*, 110–120. [CrossRef]

28. Rahani, E.K.; Kundu, T. Gaussian DPSM (G-DPSM) and Element Source Method (ESM) modifications to DPSM for ultrasonic field modeling. *Ultrasonics* **2011**, *51*, 625–631. [CrossRef] [PubMed]

29. Placko, D.; Bore, T.; Kundu, T. Distributed point source method for modeling and imaging in electrostatic and electromagnetic problems. In *Ultrasonic and Electromagnetic NDE for Structure and Material Characterization: Engineering and Biomedical Applications*; Kundu, T., Ed.; CRC Press: Boca Raton, FL, USA, 2004; pp. 249–294.

<figure><img src="image_13.png" /></figure>