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

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 planer transducer was compared with the classical DPSM formulation and analytical solution. The results show a significant improvement in the near field computation.


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][2][3][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][10][11][12][13][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][16][17][18][19][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][22][23] has been developed relatively recently for solving electrostatic, electromagnetic, 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.

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.

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. 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.

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.

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 rs.Spheres touch the front face of the transducer.Thus, the point sources are located at a distance of rs 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. From [3] the pressure and velocity fields in a fluid at point x at a distance rm from the m-th point source of strength Am can be written as: The pressure field: The velocity in the radial direction, at a distance r from the m-th point source: Therefore, the z-component of the velocity is: 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: The velocity in the radial direction, at a distance r from the m-th point source: Therefore, the z-component of the velocity is: 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: 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: 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 , . . ., 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: 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: where V S is the (M × 1) vector of the velocity components at M number of target points on the transducer surface, and A S is the (M × 1) vector containing the strengths of M number of point sources.M SS is the (M × M) square matrix relating the two vectors V S and A S : 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: The matrix M SS is obtained from Equation (7): (11) where: Equation ( 8) can be inverted to obtain the point source strengths: Instead of velocity, if pressure is specified on the transducer surface, then the source strength is obtained from the following equation: 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 where V * S is a (3M × 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 × 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: For a set of N target points the vector P T contains N pressure values at N target points, so its size is (N × 1).V T is a (N × 1) vector containing z-components of velocity at N target points, and V * T is a (3N × 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.

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)).

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)).

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)).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: The i-th family strength is denoted by λ i and its distribution is denoted by the vector q.The distribution vector q is the same for all M families but the strength λ 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: where m i jk 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: Or: Therefore: It should be noted that both Equations ( 13) and ( 20) require inversion of the M × 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 1 q 2 q 3 • • • q (F−1) q F T (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 1 11 . . .m 1 1F of Equation ( 18) tends to the same value, which corresponds to the m 11 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 1 q 2 q 3 • • • q (F−1) q F .

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 3 and P-wave speed (c f ) = 1480 m/s.To obtain the relative strength variation vector ⋯ (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 … of Equation ( 18) tends to the same value, which corresponds to the 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 ⋯ .

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  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.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 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: where V(i)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)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 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: where V(i) 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) 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 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: where V(i)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)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.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.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.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, rs/3, 2rs/3 and rs) from the transducer surface.These figures are shown in Figure 8.It should be noted that rs is the radius of the source sphere, as shown in Figure 1.Therefore, rs 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 rs, or when they are located at a distance 2rs 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 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 s /3, 2r s /3 and r s ) from the transducer surface.These figures are shown in Figure 8.It should be noted that r s is the radius of the source sphere, as shown in Figure 1.Therefore, r s 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 s , or when they are located at a distance 2r s 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.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]: 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 z ) at the transducer surface is v 0 .The wave number k f is equal to ω/c f , where c f 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.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.3rs to 1.7rs with an interval of 0.2rs 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 (vz = 1) than the classical DPSM-generated results.It should be mentioned here that Figures 5-9 were generated for σ = 1.5rs.
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 < rs) 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.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 s to 1.7r s with an interval of 0.2r s 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 z = 1) than the classical DPSM-generated results.It should be mentioned here that Figures 5-9 were generated for σ = 1.5r s .
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.
formulation are much closer to the true solution (vz = 1) than the classical DPSM-generated results.It should be mentioned here that Figures 5-9 were generated for σ = 1.5rs.
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 < rs) 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.

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 surfacebecause DPSM 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

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 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.

Figure 1 .
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.

Figure 1 .
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.
Appl.Sci.2016, 6, 302 6 of 15 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 Figure2although their average value matches with the prescribed velocity value of the transducer face.

Figure 2 .
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.

Figure 3 .
Figure 3.Quantum source clusters or source families replace the individual point sources in the new formulation for modeling the transducer.

Figure 2 .
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.
Appl.Sci.2016, 6, 302 6 of 15 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 Figure2although their average value matches with the prescribed velocity value of the transducer face.

Figure 2 .
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.

Figure 3 .
Figure 3.Quantum source clusters or source families replace the individual point sources in the new formulation for modeling the transducer.

Figure 3 .
Figure 3.Quantum source clusters or source families replace the individual point sources in the new formulation for modeling the transducer.
. 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 3 and P-wave speed (cf) = 1480 m/s.

Figure 4 .
Figure 4. Normal velocity variation (or radiation pattern) from one source family.

Figure 4 .
Figure 4. Normal velocity variation (or radiation pattern) from one source family.

Figure 5 .
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).

Figure 5 .
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).

Figure 5 .
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).

Figure 6 .
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.

Figure 6 .
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.

Figure 7 .
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.

Figure 7 .
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.

15 Figure 8 .Figure 8 .
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.Classical and family ESD DPSM results are then compared with the analytical solution in Figure9.For a circular transducer of radius of the analytical expression of the pressure field along the central axis is given by[23]:, ω ρ exp exp(22)

Figure 9 .
Figure 9. Pressure variation along the central axis in front of the transducer: dashed-dotted linetheoretical 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.

Figure 9 .
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.

Figure 10 .
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.

15 Figure 10 .
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.

Figure 11 .
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 Figure1; solid lines-computed results using source clusters, as shown in Figure3, for three different radiation patterns.

Figure 11 .
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 Figure1; solid lines-computed results using source clusters, as shown in Figure3, for three different radiation patterns.