Analysis and prediction of shock formation in acoustic energy transfer systems

Acoustic energy transfer (AET) is a transformative contactless energy transfer (CET) technology that utilizes acoustic waves to transfer energy between piezoelectric transducers. AET systems have been shown to outperform conventional electromagnetic CET technologies^1–6 in critical applications.^7,8 In particular, AET has been proposed to recharge and communicate with low-power (e.g., $1μW−−10mW$⁠) implanted medical devices,^9–11 which eliminates the need for surgery to replace batteries^11–17 and to develop battery-free underwater sensing networks to observe ocean conditions, track migration and habitats of marine animals, and monitor oil spills.^8,18–22 These applications present the need in developing mathematical and numerical models capable of assessing the efficiency of AET systems.^23–37 In general, most of the current approaches neglect nonlinear effects associated with acoustic wave propagation³⁸ and the electro-elastic response of the piezoelectric transmitters and receivers.^39,40 On the other hand, these effects become significant as the source strength is increased to enable higher energy transfer. As such, there is a need to expand the analysis capabilities to models that account for nonlinear effects and investigate their impact on the efficiency of energy transfer systems. In our previous work, we investigated theoretically the effect of material nonlinearity of a piezoelectric receiving disk on the energy transfer. We showed that the material nonlinearity can shift the optimum load resistance and that the shift is a function of the source strength.²⁹ We have also shown experimentally that the interplay of all the nonlinearities and the standing wave effects between the transmitter and receiver in an AET system can manifest themselves in a complex manner and have an impact on the energy transfer efficiency.³⁶ 

One consequence of nonlinear acoustic wave propagation is exhibited by the distortion of its waveform due to a difference in the traveling speeds of the compression and rarefaction parts of the wave. In the frequency domain, this distortion is interpreted as energy transfer from the fundamental wave frequency to its higher harmonics. The accumulation of the distortion effect as the wave travels results in a discontinuity, referred to as a shock.³⁸ The occurrence of a shock is associated with a significant loss in energy that is proportional to the cube of the difference in pressure across the discontinuity. This loss compounds as the wave propagates further, resulting in further reduction of the acoustic power.⁴¹ It is relevant to point out here that in a weakly viscous fluid such as water and air, shocks occur before the attenuation is significant.⁴² In other words, attenuation effects can be neglected, and a lossless second-order wave equation can be used to analyze the nonlinear wave propagation. In the case of a finite amplitude plane wave, the amplitudes of the higher harmonics grow at the expense of the amplitude of the fundamental or excitation frequency up to the initial shock formation location $x¯$⁠. Beyond $x¯$⁠, all components decay due to the energy dissipated as a consequence of the shock propagation.⁴² In the context of AET, the component of pressure at the excitation frequency $pω$ that is generated by the transmitting disk operating under a high excitation voltage decreases as the distance from the disk is increased due to diffraction and transfer of energy to higher harmonics. However, the total acoustic power remains conserved up to $x¯$⁠. Beyond $x¯$⁠, in addition to the transfer of energy and diffraction, the decrease in $pω$ will be compounded by additional losses in energy due to the formation and propagation of shocks. In AET systems, the power transfer efficiency depends mostly on the pressure associated with the excitation frequency, $pω$⁠, as the higher harmonics do not necessarily coincide with the higher modes of the receiver.³⁷ As such, the power transfer efficiency will be significantly compromised beyond $x¯$⁠, which renders this distance as an essential design parameter for high-intensity AET power transfer.

Analytical expressions for $x¯$ are readily available for plane waves, spherical waves, and on the axis of a focused Gaussian beam,^42,43 which is not the case in the AET system where the disk undergoes transverse deformations. Several numerical and analytical studies investigated nonlinear wave propagation and shock characteristics of acoustic waves generated by disks.^44–51 The numerical simulations provide an accurate description of the response. However, they are computationally expensive, especially when solved in the time-domain. The objective of this effort is to develop an analytical approach to predict shock formation associated with a propagating acoustic wave generated by a vibrating disk with an arbitrary transverse displacement. We consider an axisymmetric-baffled-vibrating piezoelectric disk and use the Westervelt equation to investigate the associated nonlinear wave propagation. In particular, we scale the governing equation and boundary conditions with $ϵ$ and obtain analytical expressions for the $ϵ$-order and $ϵ2$-order solutions using the Rayleigh integral. Next, we follow the work of Kelly and Nayfeh⁵² with some modifications to eliminate the secular terms by implementing the method of renormalization^53,54 and obtain a uniformly valid solution of the Westervelt equation. We validate the predictions of the analytical approach with higher fidelity finite element simulations and previously published experimental results.

An axisymmetric-baffled piezoelectric disk with thickness $h$ and radius $a$ in a semi-infinite fluid medium is considered to analyze the nonlinear wave propagation of its generated acoustic wave as shown schematically in Fig. 1. The schematic defines Cartesian $(x,y,z)$⁠, cylindrical $(rp,ψ,z)$⁠, and spherical coordinate $(r,θ,ψ)$ systems with the origin at the center of the disk $Oc$⁠. An additional spherical coordinate system $(rs,θs,ψ)$ with origin, $Os$⁠, at $z=r0$ is also defined and used in the analysis. The disk is actuated using a dynamic potential difference $V(t)$ across its flat surfaces at a frequency near that of its thickness mode. The resulting axisymmetric radial and transverse displacements are represented by $u^(rp,z,t)$⁠, and $w^(rp,z,t)$⁠, respectively. The transverse displacement of the thickness mode is chosen because it has a non-zero mean and is, therefore, favorable for generating acoustic pressure field.^34,55

The normal velocity continuity condition at the piezo-medium interface yields the necessary radiation boundary condition written as

where $w(rp,t)=w^(rp,z,t)|z=0$ and $ϕ(rp,z,t)$ is the velocity potential of the fluid. This boundary condition is an approximation because the displacement, written using the Lagrangian description, would have to be mapped to an equivalent Eulerian description to obtain the exact boundary condition. Still, this approximation is acceptable because $rp⋙u^(rp,z,t)|z=0$⁠, except at $rp=0$ where $u^(rp,z,t)|z=0=0$⁠; consequently, the effects of the approximation are insignificant. Because the modal deformation of the disk can take a complicated shape and does not possess a closed-form expression,⁵⁵ we write $w(rp,t)$ as

where $ϵ$⁠, $c0$⁠, and $ω$ are, respectively, the acoustic Mach number, velocity of sound in fluid, and excitation frequency, and $G(rp)$ and $⊘(rp)$⁠, are, respectively, the displacement amplitude and phase as a function of the distance $rp$⁠.

The lossless form of the Wesetervelt equation is used to analyze the nonlinear and non-planar wave propagation of the acoustic pressure $p(rp,z,t)$⁠. It is written as

where $ρ0$ and $β$ are, respectively, the density and the coefficient of nonlinearity of the fluid. It is relevant to note that the lossless form is valid if the analysis is restricted to distances smaller than $1/α−x¯ϵ$⁠, where $α$ is the absorption loss parameter of the fluid and $x¯ϵ$ is the shock formation distance for a plane wave.⁴² To solve Eq. (3) subjected to the boundary condition given by Eq. (1), we asymptotically expand the pressure as

and the velocity potential as

Substituting the assumed expansions of $p(rp,z,t)$ and $ϕ(rp,z,t)$ into Eqs. (1) and (3) yields the $ϵ$-order equation,

the corresponding boundary condition,

and the $ϵ2$-order equation,

From Rayleigh’s or the Kirchhoff–Helmholtz (KH) integral, the pressure generated due to a baffled-vibrating surface of surface area $S0$ and transverse velocity $U0eiωt$ is given by^56,57

where $R$ is the distance between the observation point and an infinitesimal point on the surface (⁠$AP$ in Fig. 1) and $k=ω/c0$⁠. As such, the $ϵ$-order solution is given by

The $ϵ$-order solution in Eq. (8) contains the independent variables $rp$ and $z$ inside the integral. Hence, a partial differential equation with an integral forcing function needs to be solved to determine the $ϵ2$-order solution. Such a solution is not straightforward. To simplify the analysis, we rewrite $p1(rp,z,t)$ in terms of an infinite series as suggested by Hasegawa et al.,⁵⁸ using the spherical coordinate system $(rs,θs,ψ)$ with the origin at $z=r0$ (see Fig. 1), as

where

and

Following Hasegawa et al.,⁵⁸ Eq. (9a) is modified to

where

$r12=rp2+r02$⁠, $ra2=a2+r02$⁠, $Pn$ is the Legendre function, and $hn(2)$ is the spherical Hankel function of the second kind, which is defined by $hn(2)=jn−iyn$⁠, where $jn$ and $yn$ are, respectively, the spherical Bessel functions of the first and second kind. For subsequent analysis, the $ϵ$-order solution is converted from an exponential to a trigonometric form and written as

where

Substituting the $ϵ$-order solution in the $ϵ2$-order, Eq. (6c) yields

where

and

To obtain the analytical expression for the $ϵ2$-order solution, the product of Legendre functions is rewritten as⁵² 

where the coefficient $κqnm$ is determined from the orthogonality condition of the Legendre functions as

Substituting Eq. (11) into Eq. (10a), we obtain

As detailed in the  Appendix, the solution of $ϵ2$-order equation (13), $p2(rs,θs,t)$⁠, is obtained using the separation of variables. The final solution of the pressure is then written as

The $ϵ$-order and $ϵ2$-order terms in Eq. (14) decay at different rates with respect to $rs$⁠. If left untreated, the $ϵ2$-order terms can become significant and contradict the scaling condition that $ϵ$-order terms $>>ϵ2$-order terms. To solve this contradiction and eliminate terms leading to this contradiction, we implement the method of renormalization and introduce the coordinate transformation,⁵² 

Substituting the above transformation into Eq. (14) and applying the Taylor expansion yields

From Eq. (16), the necessary condition to eliminate the secular terms is

The transformation $f$ becomes singular unless there is a particular relation between the phases of $p2$ and $∂p1/∂rs|rs=η$⁠. To examine the phase relation, we represent $p2$ and $∂p1/∂rs|rs=η$⁠, respectively, as $p2(η,θs,t)=p¯2(η,θs)cos(2ωt−2kη+⊘p2)$ and $∂p1/∂rs|r=η=p¯1′(η,θs)cos(ωt−kη+⊘p1′)$⁠, where $p¯2(η,θs)$⁠, $⊘p2(η,θs)$⁠, $p¯1′(η,θs)$⁠, and $⊘p1′(η,θs)$ are, respectively, the amplitude and phase of $p2$ and the amplitude and phase of $∂p1/∂rs|rs=η$⁠. From Eq. (17), $f(η,θs,t)$ is then determined as

It is evident from Eq. (18) that singularities arise in $f$ because of the tangent functions unless $⊘p1′−⊘p2/2=π/4$⁠. However, $⊘p1′$ and $⊘p2$ do not always hold such a relation. To eliminate this singularity, we rewrite $p2(η,θs,t)$ as

Then, eliminating the first part on the right hand side of Eq. (19) and using the remaining time-independent term as a feedback error to the $ϵ$-order solution yields

where

From Eq. (22), it is noted that the transformation $f$ is independent of $ϵ$⁠. This is a powerful consequence of the method of renormalization as once $f$ is determined from $p1$ and $p2$⁠, which are again independent of $ϵ$⁠, it can be used for any value of $ϵ$⁠.

Equations (20)–(22) constitute the solution of the nonlinear wave propagation of the acoustic pressure generated by a vibrating disk with transverse excitation according to Eq. (2). For a given $rs$⁠, $θs$⁠, and $t$⁠, one needs to first evaluate $f$ from Eq. (22) and then solve for $η$ using Eq. (21). The value of $η$ can then be used to determine the nonlinear pressure from Eq. (20). As the wave propagates in the medium, Eq. (21) will eventually yield multiple solutions for $η$ due to the cumulative nature of the nonlinearity.^38,49 The first location is where multiple solutions occur or when $∂p/∂rs=∞$ is the shock location.^{38,52,59–61} Beyond this location, the solution is not valid, and a shock fitting criteria such as an equal-area rule should be used.^38,53

The efficacy of the analysis presented in Sec. II to predict the nonlinear wave propagation and shock formation is assessed by comparing its predictions with those from (a) higher fidelity Finite Element (FE) simulations and (b) previously published experimental results.

A piezoelectric disk with radius $a=5$ mm and thickness $h=2$ mm made of PZT-5H material and submerged in water was considered for validation using numerical simulations. The material and acoustic properties of the disk and fluid are presented in Table I. We performed a series of axisymmetric frequency and time-domain FE simulations in COMSOL Multiphysics to obtain the acoustic radiation characteristics of the baffled-disk configuration. A quarter-circular fluid domain of radius $rf$ was considered in the simulations. A spherical wave radiation boundary condition at the outer boundary was used to simulate an infinitely propagating wave, and triangular mesh elements were used in all the domains.

In both analysis and simulations, the material and geometric nonlinearities of the disk are neglected. Furthermore, because of the presence of a baffle, only the transverse displacement of the top surface of the disk is required to determine its acoustic radiation. In the simulations, we extracted this displacement from a linear FE simulation and use it as a boundary condition in the nonlinear FE simulation for computational efficiency, realized by eliminating the multiphysics coupling between the piezoelectric and fluid domains.

Considering a maximum mesh size of $λf0/6$ in the fluid domain and 20 elements each along the radius and thickness of the disk, we solved the linear problem using the frequency domain FE solver and determined that the thickness extensional mode resonant frequency of the disk is $f0=1.005$ MHz (see p. 108 in Guo⁶²), where $λf0=c0/f0$⁠. The corresponding Rayleigh far-field distance is $rfar=12ka2=53.3$ mm. The response and radiation characteristics of the disk for an excitation voltage amplitude corresponding to $ϵ=10−5$ and an excitation frequency of $f0$ are presented in Fig. 2. In particular, Figs. 2(a)–2(d) show, respectively, the sound pressure level (⁠$pref=1$ Pa) in the fluid domain, deformation of the disk, transverse displacement of the top surface of the disk, and the beam pattern at different radial distances from the origin [see schematic in Fig. 2(a)]. The abnormal higher deformation amplitude on the bottom surface compared to the top surface of the disk in Fig. 2(b) is due to the absence of radiation damping unlike the top surface. From Fig. 2(d), we note that except for $r=0.125rfar$⁠, the pressure along the $z$-axis is maximum for any $r$ value. Consequently, we expect the shock to occur first on the $z$-axis. We also note that $r=0.125rfar$ is a local minimum similar to those observed in the near-field of axial pressure generated by a piston.^46,57,63 Because we expect the shock to occur along the $z$-axis first, we limit further analysis only to the axial pressure distribution.

Next, the transverse displacement presented in Fig. 2(c) is used to obtain axial $p1$ and $p2$ according to the solution in Sec. II. Alternatively, $p1$ and $p2$ could be determined by solving the $ϵ$-order and $ϵ2$-order equations and the $ϵ$-order boundary condition in a frequency domain FE solver. A pressure release boundary condition on the disk should be used to solve the $ϵ2$-order equations to ensure that $p2$ vanishes on the surface of the disk. We note that $p2$ evaluated using the two approaches should differ only by the non-secular terms (NSTs) that were neglected in the analysis. The axial $p1$ and $p2$ evaluated from the model for $n=50,r0=−50/k$⁠, labeled as Series $n=50$⁠, and by the frequency domain FE simulations are presented, respectively, in Figs. 3(a) and 3(b). A mesh size of $λf0/12$ was used in the FE simulations. Figure 3(a) shows an excellent agreement between $p1$ evaluated by the analysis and numerical simulations except at close distances to the disk. The disagreement in this region is due to limiting the series and can be minimized by choosing higher values of $n$ and $r0$⁠. From Fig. 3(b), we note a close agreement in the amplitude and a small phase shift. We attribute this shift to errors induced by choosing a finite value of $n$ and to the contribution of the NST. It is relevant to point out that $p1$ and $p2$ determined using FE are more accurate and that eliminating $p2$ that includes NST in the method of the renormalization will result in a transformation that will account for the complete particular solution of the $ϵ2$-order equation. As for computational cost, the evaluation of $p1$ and $p2$ from the analysis requires the analytical expressions of integrals presented in the  Appendix and storing, reusing, and manipulation of these expressions, which can become computationally expensive, especially for higher values of $n$⁠. On the other hand, the linear frequency domain FE simulations to determine $p1$ and $p2$ is more efficient as it involves a matrix inversion. Moreover, the FE solver can be easily extended to complicated geometries.

Having determined $p1$ and $p2$⁠, we compare next the nonlinear wave propagation predicted using Eqs. (20)–(22) with that determined from a nonlinear time-domain FE simulation for three different excitation levels, namely, $ϵ=5×10−3$⁠, $ϵ=2.5×10−3$⁠, and $ϵ=1.25×10−3$⁠. Toward that objective, we constructed a fluid domain of $rf=2.25×x¯ϵ$ that is excited by the transverse displacement presented in Fig. 2(c). A maximum mesh size of $λf0/36$ was selected so that the response at $6f0$ can be captured. For the sake of completeness, the Westervelt equation with dissipation was simulated in time-domain using COMSOL Multiphysics by choosing the diffusivity of the sound as $δ=1.34×10−6m2/s$⁠. However, the absorption is not expected to have a significant impact on the response because the characteristic absorption distance [$1/α=2c03/(δω2)$] for $f0$ is about 120 m, which is significantly larger than $2.25×x¯1.25×10−3(≈0.1m$⁠). The acoustics (Westervelt) model along with the shock-capturing stabilization scheme with $q=1.35$ and $κ=0.00015$ was used to solve for the nonlinear wave propagation in this study.

The steady-state response characteristics for $ϵ=5×10−3$ obtained using time-domain FE simulation and from Eqs. (19)–(22) using $p1$ and $p2$ as determined from the analysis (perturbation and series, $n=50$⁠) and frequency domain FE simulations (perturbation and freq. FE) are compared next. The plots presented in Figs. 4(a)–4(d) show, respectively, the steady-state axial pressure waveform at time $t=1/(2f0)$⁠, the mean intensity, the amplitude of the pressure component with frequency $ω=2πf0$⁠, and the amplitude of the pressure component with frequency $2ω$⁠. Figures 4(b) and 4(c) also show, respectively, the mean intensity and pressure when $β=0$⁠, i.e., for the linear case. Based on the closeness of the responses predicted when $p1$ and $p2$ are determined from the analysis and from the FE simulations, it can be concluded that the contribution of NST is not significant. It can also be concluded that the transformation does an excellent job in predicting the pressure distributions at closer distances and that the discrepancies grow as the distance increases. The discrepancy is because the $ϵ3$-order secular terms become dominant as the distance increases to the point where one needs to carry the method of renormalization to $ϵ3$-order to maintain the accuracy. Figures 4(a)–4(d) also show regions where the transformation yields multiple solutions, which invalidates the solution, and as such, no solution is shown. As noted in previous studies,^38,53 these regions correspond to the shock formation, and a shock fitting criteria should be used to obtain the true waveform in these regions. From the results, the transformation yields multiple solutions in a small region around $r/x¯5×10−3=0.5$ and when $r/x¯5×10−3>1.13$⁠. However, from Fig. 4(b), the mean intensity predicted by the time-domain FE simulation deviates from the linear case only beyond $r/x¯5×10−3=1.13$⁠. This deviation is a consequence of loss of energy due to the formation and propagation of shock.^38,46 As such, we hypothesize that $r/x¯5×10−3=1.13$ is the shock formation location for $ϵ=5×10−3$⁠. As expected, an accelerated decrease in $pomega$ is seen after the shock formation at $r/x¯5×10−3=1.13$ in Fig. 4(c). To confirm that a shock did not occur around $r/x¯5×10−3=0.5$⁠, we present the time series and Fourier coefficients of the pressure at $r/x¯5×10−3=0.48$ as obtained from the time-domain FE simulation, respectively, in Figs. 5(a) and 5(b). The Fourier component of pressure at $2f0$ is comparable with that at $f0$ in Fig. 5(b), and the corresponding time series presented in Fig. 5(a) does not indicate a shock, confirming that no shock occurs in the region around $r/x¯5×10−3=0.5$⁠. We attribute this anomaly predicted by the model to the local minima in the amplitude of pressure, which is essentially a singularity. As such, multiple solutions are also associated with local minima of the pressure amplitude in the near field. To the knowledge of the authors, this conclusion has not been made in previous investigations as the pressure fields analyzed in those studies did not exhibit local minima.

The steady-state axial pressure waveform at time $t=0$⁠, mean intensity, the component of pressure at $ω$⁠, and the component of pressure at $2ω$ akin to Figs. 4(a)–4(d) for the case of $ϵ=2.5×10−3$ are, respectively, presented in Figs. 6(a)–6(d). The agreement of the results obtained from the analysis and the time-domain FE simulation is similar to that observed in Figs. 4(a)–4(d). The only major difference is the absence of an anomaly, which suggests that the occurrence of anomaly also depends on the magnitude of $ϵ$⁠. More importantly, the results predicted by the analysis suggest that the shock formation location on the $z$-axis is $r/x¯2.5×10−3=1.04$⁠, which is also the point after which the mean intensity predicted by the time-domain FE simulation deviates from the linear case. The characteristic shock formation distance for $ϵ=2.5×10−3$ is $x¯2.5×10−3=26.07$ mm.

The results for the case of excitation amplitude $ϵ=1.25×10−3$ are presented in Figs. 7(a)–7(d). Inspecting the steady-state axial pressure waveform at time $t=0$ predicted by the time-domain FE simulation in Fig. 7(a), we note that a shock does not take place in the range of $0<r<2.25×x¯1.25×10−3$⁠, where $x¯1.25×10−3=52.13$ mm. This result can also be inferred from Fig. 7(b), which does not show any deviation of the mean intensity predicted by the time-domain FE simulation from the linear case. Consequently, the deviation of $pω$ obtained from the time-domain FE simulation from the linear case in Fig. 7(c) is very small. The results predicted by the analysis suggest that the shock formation location on the $z$-axis is around $r/x¯1.25×10−3=1.64$⁠, which falls in the far-field region of the disk. This inconsistency is again a consequence of implementing the method of renormalization only up to the O(⁠$ϵ2$⁠) order, as discussed earlier.

The efficacy of the developed model in predicting the nonlinear wave propagation and the shock formation location is also validated by comparing its predictions with the experimental data of Nachef et al.⁶³ Nachef et al. considered a piezoelectric ceramic disk having a radius of 25 mm and measured the pressure waveforms over axial distances between 50 mm and 400 mm and lateral distances between 0 mm and 50 mm for different source strengths at an excitation frequency of $f0=1$ MHz. These waveforms were later treated as a benchmark by Khokhlova et al.⁴⁶ to validate their numerical models for $p0=0.2,0.64,$ and $1.43$ MPa where $p0=ϵρ0c02$ is the amplitude of the pressure at the source. They identified that the shock formation distances on the axis are, respectively, $0.1rfar$ and $0.28rfar$ for $p0=1.43$ MPa and $p0=0.64$ MPa where the Rayleigh far-field distance $rfar=1156.63$ mm. They concluded that a shock did not occur on the axis before 400 mm⁴⁶ when $p0=0.2$ MPa. The relevant material and geometric properties were identified as $ρ0=1000kg/m3$⁠, $c0=1500$ m/s, and $β=3.5$ and the equivalent aperture of the circular piston as $a=23.5$ mm.⁴⁶ 

From the discussion in Sec. III A 2, $p1$ and $p2$ can be evaluated either by the analysis or from appropriate frequency domain linear FE simulations. For the sake of simplicity, we evaluated $p1$ and $p2$ using the FE simulations to determine the nonlinear wave propagation generated by a baffled piston. Figures 8(a) and 8(b) show, respectively, the comparisons of Fourier components of pressure at $ω$ and $2ω$ predicted using Eqs. (20)–(22) with the linear case (⁠$β=0$⁠) and corresponding experimental data when $p0=0.2$ MPa. The corresponding results for $p0=0.64$ MPa and $p0=1.43$ MPa are presented, respectively, in Figs. 8(c) and 8(d) and Figs. 8(e) and 8(f). From the results, we note that the shock formation distances predicted by the method of renormalization are $r/rfar=0.25$ and $r/rfar=0.1$⁠, respectively, for $p0=0.64$ MPa and $p0=1.43$ MPa and that a shock did not occur in the range considered for $p0=0.2$ MPa. These findings are very close to the shock formation distances identified by Khokhlova et al.,⁴⁶ which validate the analysis.

Various effects determine the desired distance in AET systems that delivers the optimum power transfer. The output power fluctuates as the separation distance between the transmitter and receiver is increased due to the spatial acoustic resonances. These fluctuations are strong at closer distances and often become weak with an increase in the separation distance as the influence of the reflected waves diminishes with an increase in the separation distance. When these fluctuations are strong, constructive interference positions become the desired distances for AET. However, for a finite-sized transmitter and receiver disks, the general trend in the output power decreases with an increase in the separation distance. On the other hand, if the transmitter’s excitation amplitude is strong enough to trigger the wave propagation nonlinearities, we note that the formation and propagation of shock are associated with a loss in energy, thereby compromising the output power transfer beyond the initial shock formation distance. As such, the knowledge of the shock formation distance is crucial for AET systems.

The shock formation distance of a piezoelectric transmitter disk with arbitrary deformation depends on the excitation amplitude, the operation frequency, deformation of the disk, the disk’s diameter, and the material properties of the acoustic medium such as density, speed of sound, and the coefficient of nonlinearity. A closed-form expression that estimates the shock formation distance is not available in the literature. The analysis developed in this work provides a quick estimate that can aid the design and operation of efficient AET systems. The analysis is based on the asymptotic expansion of the governing equation. The method of renormalization is used to predict the nonlinear wave generated by a finite amplitude baffled circular disk for a specified deformation profile. The lossless form of the Westervelt equation and the normal velocity continuity boundary condition were scaled with a non-dimensional parameter related to the amplitude of the transverse displacement of the disk. They were expanded and solved by introducing a coordinate transformation to eliminate the secular terms, which resulted in singularities that required redefining the relationship between the pressure components. We validated the analysis by comparing its predictions of the shock formation distance under different excitation levels with higher fidelity nonlinear finite element simulations and previously published experimental results. We also demonstrated the versatility of the analysis by showing that the $ϵ$-order and $ϵ2$-order solutions can be determined either from analytical expressions or from linear frequency domain finite element simulations. We showed that the analysis accurately predicted the nonlinear wave propagation and shock formation distance in the near-field of the disk and argued that the accuracy level could be increased by extending the analysis to $ϵ3$ or higher orders. We found out that the transformation yields multiple solutions at the shock location and at local minima that occur in the near-field of the pressure. We also showed an accelerated decrease in the power of the excitation frequency beyond the shock formation. Of particular importance is the validity of the solution for any excitation level because the $ϵ$ and $ϵ2$ order solutions are independent of that level.

This work was supported by the National Science Foundation (NSF) (Grant No. ECCS-1711139), which is gratefully acknowledged. The authors would also like to thank Professor Saad Ragab (Virginia Tech) for providing valuable feedback in developing the model.

The authors declare that they have no conflict of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To determine the particular solution of Eq. (13), we use the separation of variable technique and assume the solution as

Substitution of Eq. (A1) in Eq. (13) yields

where

The homogeneous solutions of Eqs. (A2a) and (A2b) are $X1(rs)=Ψ1(rs)=jq(2krs)$ and $X2(rs)=Ψ2(rs)=yq(2krs)$⁠, which are now used to determine its particular solution by using the variation of the parameter technique as