© C. Qi et al., Published by EDP Sciences, 2021

1 Introduction

The active control of acoustic fields has been extensively explored in the past decades and is an emerging research area in modern acoustics. Compared with passive control schemes, active control techniques are more suited for low frequencies, as the barriers used in passive control schemes are relatively large and ineffective at low frequencies [1]. Instead of using the interaction between sound and specific materials to control the acoustic field, active control strategies characterize a source so that it is capable of approximating given field patterns in prescribed exterior regions. The current literature has significantly addressed the idea of the active control of Helmholtz scalar fields in broad applications. These include, but are not limited to, active noise control [2–4], personal sound zones or multizone sound reproduction [5–7], active acoustic cloaking [2, 8–10], remote sensing [11, 12] and metamaterial design [13–15]. Active sound control techniques are becoming increasingly ubiquitous to enhance sound-based systems [16–21]. Forward and inverse problems for sound in underwater environments have been widely studied in the literature (see monographs [22–24] and references therein). Several comprehensive reviews are available in [1, 25, 26], discussing the technicalities and recent advances in various active control schemes.

Directional far-field control can potentially play an important role in some applications with far-field manipulation involved. In this context, in [27, 28], the authors investigated the directional acoustic manipulation via multi-phase forehead structure of the porpoises to realize beamforming. In the same paradigm, the works [29–31] made use of loudspeaker arrays to achieve a directional source.

The majority of the active control strategies have been focused on the method of pressure matching (PM) [25]. The PM approach aims to match the target field pattern in the given region with minimum error. In this approach, the active control problem is cast as an inverse source problem (ISP).

Unlike in free space, the active sound control in underwater environments is more challenging. It makes use of additional boundary conditions and the approximated Green’s function, which adds a layer of complexity than free space scenario. We mention here the works [32, 33] where the authors develop a single-mode excitation with a feedback control algorithm to realize both near and far-field sound control.

A general ocean environment can be modeled as a horizontally stratified waveguide [23, 24, 34]. In general, it is fairly difficult to find the analytical fundamental solution and thus, in this paper we follow the paradigm proposed in [22] and consider a simpler marine environment modeled as a shallow water or a homogeneous finite-depth ocean. Our sensitivity analysis builds up on the numerical framework developed in [20, 21, 35]. We use the associated Green’s function to represent the solution of the Helmholtz equation and employ the integral equation (IE) method to formulate the forward propagator. The method of moments (MoM) approach [36] is used to reduce the original integral equation to a discrete linear system. Then a Tikhonov regularization scheme with the Morozov discrepancy principle is applied to solve the resulting system of equations. In the underwater environment, additional boundary conditions, including pressure release boundary on air-water interface and sound hard boundary on the sea floor, were imposed. Consequently, the Green’s function was modified in the formulation of the forward propagator. We used the normal mode representation to formulate the Green’s function in the homogeneous finite-depth ocean [22, 23].

In this paper, we present a detailed sensitivity study for the problem of controlling three-dimensional scalar Helmholtz fields in several prescribed exterior regions while maintaining desired far-field pattern values in given fixed directions. Moreover, we also look at the process error, i.e., small noise that may be present in the feeding current when the source is physically implemented. We discuss the feasibility of the active scheme (in terms of power budget, control accuracy and process error) with respect to variations in frequency, the distance between the control region and the source, the mutual distance between the control regions, and the control region size. By a feasible source, we mean a source that requires low power (within the dB levels of common sound sources such as loudspeakers), achieves the control effects with high accuracy (errors of less than 10%) and with high tolerance to feeding noise.

The rest of this paper is organized as follows. In Section 2, we formally describe the problem and provide relevant theoretical results obtained in [35]. Section 3 shows the numerical results and sensitivity analysis in free space. In Section 4, the numerical results and sensitivity analysis in the shallow water environment are presented. Finally, we conclude the paper with some remarks in Section 5.

2 Theory

2.1 Problem formulation

In this section, we present a general description of the active manipulation scheme for Helmholtz fields proposed in our previous works. The unified functional and numerical framework have already been discussed in [20, 21, 35]. We shall briefly recall several essential theoretical results and describe some geometric configurations of interest.

The problem is to characterize a source (modeled as surface pressure or surface normal velocity) so that its generated field approximate some prescribed fields in several exterior regions of interest while maintaining desired patterns in several given far field directions. In this paper, the active field manipulation scheme is explored in both free space and homogeneous ocean with a constant depth. The problem geometry in these two environments are sketched in Figures 1 and 2. Although the theoretical discussion in [20, 21, 35] indicates that an arbitrary number of source regions, exterior control regions, and far-field directions can be considered in the active control scheme, here we shall only consider a single source D_a, two control regions D₁, D₂, and two far-field directions x₁, x₂ for illustrative purposes. A single source D_a ⋐ ${\mathbb{R}}^3$ is modeled as a compact region in both free space and homogeneous ocean. The control regions D₁ and D₂ are mutually disjoint smooth domains, i.e., D₁ ∩ D₂ = ∅. We also assume that the control regions are well-separated from the source region, i.e., (D₁ ∪ D₂) ∩ D_a = ∅. Furthermore, we consider two distinct directions x₁ and x₂ representing the far-field directions of interest, which can be arbitrary.

Figure 1 Sketch of the problem geometry showing the near controls D1, D2 and the far field directions x1 and x2 in free space.

Figure 2 Sketch of the problem geometry showing the near controls D1, D2 and the far field directions x1 and x2 in homogeneous ocean with constant depth.

Mathematically the problem is to find the boundary input on the source, either a Dirichlet input data p (pressure) or a Neumann input data v_n (normal velocity) such that for any desired field f = (f₁, f₂) on the control regions D₁, D₂ and prescribed far field pattern values f_∞ = (f_∞,1, f_∞,2), the solution u of the following exterior Helmholtz problem:

$$\{\begin{array}{llll}{\nabla }^{2}u+k^{2}u=0 \mathrm{in} {\mathbb{R}}^3\backslash D_a,& & & \\ \nabla u·\mathbf{n}=v_n,(\mathrm{or} u=p) \mathrm{on} \mathrm{\partial }{D}_a,& & & \\ \mathrm{corresponding}\hspace{1em}\mathrm{boundary}\hspace{1em}\mathrm{conditions},& & & \\ \mathrm{suitable}\hspace{1em}\mathrm{radiation}\hspace{1em}\mathrm{condition},& & & \end{array}$$(1)satisfies the control constraints,

$$\{\begin{array}{c}{\Vert u-f_l\Vert }_{c^2\left(D_l\right)}\le \mu \hspace{1em}\mathrm{for}\hspace{1em}l=\overline{\mathrm{1,2}},\\ \left|u_{\infty }\left({\widehat{\mathbf{x}}}_j\right)-f_{\infty ,j}\right|\le \mu \hspace{1em}\mathrm{for}\hspace{1em}j=\overline{\mathrm{1,2}},\end{array}$$(2)

where 0 < μ ≪ 1 is the desired control accuracy threshold and u_∞ denotes the the far field pattern of u. In (1) and (2), n denotes the outward-pointing normal vector to ∂D_a and $\widehat{\mathbf{x}}=\frac{\mathbf{x}}{\left|\mathbf{x}\right|}$ is the unit vector along the direction x. The subscript C² denotes that the norm is defined in the space of smooth functions with continuous derivatives up to the second order and is computed as a sum of the L² norms of all the partial derivatives up to the second order. Recall that the L² norm on the space of square integrable functions on a given domain D is defined as ${\Vert f\Vert }_{L^2\left(D\right)}={\left({\int }_D\mathrm{}\left|f\right(\mathbf{x}){|}^{2}d\mathbf{x}\right)}^{1/2}$. The boundary and radiation conditions in (1) depend on the medium (environment) and will be given in Sections 2.2 and 2.3.

In [20, 21], it was demonstrated that problem (1) together with (2) admits a solution if the wavenumber k is not a resonance, i.e., −k² is neither a Neumann eigenvalue for the Laplacian in the source region nor a Dirichlet eigenvalue in the control regions. In order to ease the analysis and integral computations, our scheme makes use of a “fictitious source”, i.e., an arbitrary sphere $D_a^{\text{'}}$ compactly embedded in the actual source region D_a. In general, the physical source D_a can have any shape as long as it has a Lipschitz boundary, compactly includes the fictitious source $D_a^{\text{'}}$ and is well separated from the control regions. Meanwhile, our scheme uses slightly larger mutually disjoint regions W₁ and W₂ such that D₁ ⋐ W₁, D₂ ⋐ W₂, $W_1\cap W_2=\mathrm{\varnothing }$ and $(W_1\cup W_2)\cap D_a=\mathrm{\varnothing }$ because, as shown in [20], an accurate control in the sense of the L²-norm on ∂W₁ and ∂W₂ implies, via regularity and uniqueness results for the solution of the interior Helmholtz equation, the smooth control required in (2). As pointed out in [21, 35], within the framework mentioned above, the boundary input data, either normal velocity v_n or pressure p on the surface of the source can be characterized using a smooth density function $w\in L^2\left(\mathrm{\partial }{D}_a^{\text{'}}\right)$ such that,

$$v_n\left(\mathbf{x}\right)=\frac{-i}{\mathrm{pck}}\frac{\partial }{\partial \mathbf{n}}{\int }_{\partial D_a^{\text{'}}}^{}w\left(\mathbf{y}\right)\varphi \left(\mathbf{x},\mathbf{y}\right)\mathrm{d}{S}_{\mathbf{y}}\hspace{1em}\mathrm{and},$$(3)

$$p\left(\mathbf{x}\right)={\int }_{\partial D_a^{\text{'}}}^{}w\left(\mathbf{y}\right)\varphi \left(\mathbf{x},\mathbf{y}\right)\mathrm{d}{S}_{\mathbf{y}},$$(4)

for x ∈ ∂D_a and where ρ is the density of the surrounding environment, c is the speed of sound in the given medium and ϕ(x, y) is the fundamental solution of the 3D Helmholtz equation. In general, the solution of problem (1) can be represented by a linear combination of a single and a double layer potential illustrated in [37, 38]. For simplicity of the computations, we only use the single-layer potential operator throughout this paper.

2.2 Free space environment

In this section, the active manipulation of Helmholtz fields in a free-space environment is investigated. In this medium, problem (1) reads,

$$\{\begin{array}{c}{\nabla }^{2}u+k^{2}u=0 \mathrm{in} {\mathbb{R}}^3\backslash D_a,\\ \nabla u·\mathbf{n}=v_n,(\mathrm{or} u=p) \mathrm{on} \mathrm{\partial }{D}_a,\\ \langle \widehat{\mathbf{x}},\nabla u\left(\mathbf{x}\right)\rangle -\mathrm{iku}\left(\mathbf{x}\right)=o\left(\frac{1}{\left|\mathbf{x}\right|}\right), \hspace{1em}\mathrm{as} \left|\mathbf{x}\right|\to \mathrm{\infty }.\end{array}$$(5)

The fundamental solution in this case is given by $\varphi (\mathbf{x},\mathbf{y})=\frac{e^{\mathrm{ik}|\mathbf{x}-\mathbf{y}|}}{4\pi |\mathbf{x}-\mathbf{y}|}$. Then, the control problem is to characterize p or v_n such that the solution u of (5) satisfies (2). In free space, the solution u can be written as the image of a density function w under a forward propagator operator $\mathcal{D}$. To define $\mathcal{D}$, we make use of the near field operator ${\mathcal{K}}_l$ on the control region ∂W_l defined as,

$${\mathcal{K}}_{l}w\left({\mathbf{z}}_l\right)={\int }_{\partial D_a^{\text{'}}}^{}w\left(\mathbf{y}\right)\varphi \left({\mathbf{z}}_l,\mathbf{y}\right)\mathrm{d}{S}_{\mathbf{y}},$$(6)

where for each $l=\overline{\mathrm{1,2}}$, z_l ∈ ∂W_l and $\mathbf{y}\in \mathrm{\partial }{D}_a^{\text{'}}$. Moreover, following the derivation in [38], the far field pattern operator ${\mathcal{K}}_{\mathrm{\infty },j}$ can be defined as,

$${\mathcal{K}}_{\mathrm{\infty },j}w\left({\widehat{\mathbf{x}}}_j\right)=1\backslash 4\pi {\int }_{\mathrm{\partial }}\mathrm{}w\left(\mathbf{y}\right)e^{-\mathrm{ik}{\widehat{\mathbf{x}}}_j·\mathbf{y}}\mathrm{d}{S}_{\mathbf{y}},$$(7)

where ${\widehat{\mathbf{x}}}_j$, $j=\overline{\mathrm{1,2}}$ is the unit vector pointing in the far-field direction of interest. Hence, the overall propagator operator $\mathcal{D}$ in free space is defined as,

$$\begin{array}{cc}\mathcal{D}w\left({\mathbf{z}}_1,{\mathbf{z}}_2,{\widehat{\mathbf{x}}}_1,{\widehat{\mathbf{x}}}_2\right)=& \left({\mathcal{K}}_{1}w\right({\mathbf{z}}_1),{\mathcal{K}}_{2}w({\mathbf{z}}_2),\\ & {\mathcal{K}}_{\mathrm{\infty },1}w\left({\widehat{\mathbf{x}}}_1\right),{\mathcal{K}}_{\mathrm{\infty },2}w\left({\widehat{\mathbf{x}}}_2\right)).\end{array}$$(8)

2.3 Homogeneous ocean environment

Compared with the free space regime, the active control scheme in the homogeneous ocean with a constant depth is much more complicated. As shown in Figure 2, two additional boundary conditions must be specified. More explicitly, problem (1) now reads,

$$\{\begin{array}{c}{\nabla }^{2}u+k^{2}u=0 \mathrm{in} {\mathbb{R}}^3\backslash D_a,\\ \nabla u·\mathbf{n}=v_n,(\mathrm{or} u=p) \mathrm{on} \mathrm{\partial }{D}_a,\\ \begin{array}{c}u=0\mathrm{ }\mathrm{at}\mathrm{ }\mathrm{the}\mathrm{ }\mathrm{ocean}\mathrm{ }\mathrm{surface}\mathrm{ }z=0,\\ \frac{\mathrm{\partial }u}{\mathrm{\partial }z}=0\mathrm{ }\mathrm{at}\mathrm{ }\mathrm{the}\mathrm{ }\mathrm{ocean}\mathrm{ }\mathrm{floor}\mathrm{ }z=h,\\ \underset{r\to \mathrm{\infty }}{\lim }{r}^{1/2}\left(\frac{\mathrm{\partial }{u}_p}{\mathrm{\partial }r}-\mathrm{ik}{a}_p{u}_p\right)=0,\hspace{1em}\mathrm{for}\mathrm{ }\theta \in \left[\mathrm{0,2}\pi \right).\end{array}\end{array}$$(9)

where following the framework proposed in [22], we employ cylindrical coordinates in our analysis and the functions u_p’s represent the normal modes in the representation of u. The main control problem is then to characterize v_n or p such that the solution u of (9) satisfies (2). We employ the following Green’s representation for u,

$$u\left(\mathbf{x}\right)={\int }_{\partial D_a^{\text{'}}}^{}w\left(\mathbf{y}\right)G\left(\mathbf{x},\mathbf{y}\right)\mathrm{d}{S}_{\mathbf{y}},$$(10)

where w is the density function defined on the fictitious surface source $\mathrm{\partial }{D}_a^{\text{'}}$ and G is the associated Green’s function in the medium. For any observation point x = (r, θ, z) = (ξ, z) and source point y = (r′, θ′, z′) = (ξ′, z′), the Green’s function has the following normal mode representation [24]:

$$G(\mathbf{x},\mathbf{y})=\frac{i}{2h}\sum _{p=0}^{+\mathrm{\infty }}\mathrm{}{\varphi }_p\left(z\right){\varphi }_p\left(z\mathrm{\prime }\right)H_0^{\left(1\right)}\left(k{a}_p\right|\mathit{\xi }-\mathit{\xi }\mathrm{\prime }\left|\right),$$(11)

where $H_0^{\left(1\right)}$ is the zero order Hankel function of the first kind, ϕ_p is the pth modal solution with associated eigenvalue a_p [22–24]. These eigenvalues are,

$$a_p=\sqrt{1-\frac{(2p+1{)}^2{\pi }^2}{4k^2{h}^2}},$$(12)

while the separated modal solutions ϕ_p are given by,

$${\varphi }_p\left(z\right)=\sin \left[k\sqrt{1-a_p^2}z\right].$$(13)

In the far-field region, the field u has an asymptotic form given by [22],

$$u\left(\mathbf{x}\right)=\sum _{p=0}^N\mathrm{}\frac{1}{\sqrt{k{a}_{p}r}}{e}^{\mathrm{ik}{a}_{p}r}{g}_p\left(\theta ,z\right)+\mathcal{O}\left(\frac{1}{r^{\frac{3}{2}}}\right),\hspace{1em}\mathrm{as} r\to +\mathrm{\infty },$$(14)

where N is the number of propagating modes (i.e., the largest integer so that $a_p\in \mathbb{R}$), g_p is given by,

$$g_p(\theta ,z)=\sqrt{\frac{2}{\pi }}{\int }_{\mathrm{\partial }}\mathrm{}w\left(y\right)\cdot \left(\sum _{q=0}^{\mathrm{\infty }}\mathrm{}{e}^{-i(q+\frac{1}{2})\frac{\pi }{2}}{\alpha }_{\mathrm{qp}}(z,\theta ,r\mathrm{\prime },z\mathrm{\prime },\theta \mathrm{\prime })\right)\mathrm{d}{S}_y,$$(15)

and,

$${\alpha }_{\mathrm{qp}}(z,\theta ,r\mathrm{\prime },z\mathrm{\prime },\theta \mathrm{\prime })=\frac{i{ϵ}_q}{2h}{\varphi }_p\left(z\right)\cdot \left[\cos \left(\mathrm{q\theta }\right){\beta }_{\mathrm{qp}}\left(\mathbf{y}\right)+\sin \left(\mathrm{q\theta }\right){\gamma }_{\mathrm{qp}}\left(\mathbf{y}\right)\right],$$(16)

where ϵ₀ = 1, and ϵ_q = 2 for q ≥ 1 with,

$${\beta }_{\mathrm{qp}}\left(\mathbf{y}\right)=J_q\left(k{a}_{p}r\mathrm{\prime }\right){\varphi }_p\left(z\mathrm{\prime }\right)\cos \left(\mathrm{q\theta }\mathrm{\prime }\right),$$(17)

$${\gamma }_{\mathrm{qp}}\left(\mathbf{y}\right)=J_q\left(k{a}_{p}r\mathrm{\prime }\right){\varphi }_p\left(z\mathrm{\prime }\right)\sin \left(\mathrm{q\theta }\mathrm{\prime }\right),$$(18)

where J_p(x) is the Bessel function of the first kind of order p. Therefore, the far-field pattern in a given direction x = (1, θ, z) can be defined as [22],

$$u_{\mathrm{\infty }}\left(\widehat{\mathbf{x}}\right)=\sum _{p=0}^N\mathrm{}{g}_p(\theta ,z).$$(19)

Remark 2.1

In (5) and (9), we assume that the free space and homogeneous environment are idealized. In free space, there is no additional boundary condition considered. In the shallow water environment, only air-water surface and the sea floor are applied as two boundary conditions. These boundary conditions are used to demonstrate the sensitivity study. Adaptive or random boundary conditions and the corresponding Green’s function have to be imposed when a practical environment is considered.

Similar to the free space regime, we define a propagator operator $\mathcal{D}$ that calculates the generated field on the control regions and the far-field pattern in the given directions. For each $l=\overline{\mathrm{1,2}}$ and $j=\overline{\mathrm{1,2}}$, define,

$${\mathcal{K}}_{l}w\left({\mathbf{z}}_l\right)={\int }_{\partial D_a^{\text{'}}}^{}w\left(\mathbf{y}\right)G({\mathbf{z}}_l,\mathbf{y})\mathrm{d}{S}_{\mathbf{y}}\hspace{1em}\mathrm{and},$$(20)

$${\mathcal{K}}_{\mathrm{\infty },j}w\left({\widehat{\mathbf{x}}}_j\right)=\sum _{p=0}^N\mathrm{}{g}_p({\theta }_j,z_j).$$(21)

Then the overall propagator operator $\mathcal{D}$ is given by,

$$\begin{array}{ll}\mathcal{D}w({\mathbf{z}}_1,{\mathbf{z}}_2,{\widehat{\mathbf{x}}}_1,{\widehat{\mathbf{x}}}_2)=& \left({\mathcal{K}}_{1}w\right({\mathbf{z}}_1),{\mathcal{K}}_{2}w({\mathbf{z}}_2),\\ & {\mathcal{K}}_{\mathrm{\infty },1}w\left({\widehat{\mathbf{x}}}_1\right),{\mathcal{K}}_{\mathrm{\infty },2}w\left({\widehat{\mathbf{x}}}_2\right)).\end{array}$$(22)

2.4 Optimization scheme

In Sections 2.2 and 2.3, we already defined the propagator operator $\mathcal{D}$ that evaluates the field in the exterior control regions and the far-field pattern values. The problem (1), (2), formulated in the respective context of free space model (5) or constant-depth homogeneous ocean model (9) can be summarized as,

$$\mathcal{D}w\approx f.$$(23)

where the operator $\mathcal{D}$ is defined in (8) and (12), w is the unknown density function defined on the fictitious source and f denotes the prescribed fields in the control regions. Following the approach in [17] the density function w in (23) is determined using the method of moments (MoM) by discretizing the control regions into a discrete mesh of collocation points and w being expressed as a linear combination (with unknown coefficients) of local basis functions spanning the space of square integrable functions on $\mathrm{\partial }{D}_a^{\text{'}}$. Thus, the integral equation in (23) is reduced to a linear system,

$$\mathbf{A}{\mathbf{w}}_d=\mathbf{b},$$(24)

where w_d represents the vector of unknown coefficients in the local basis representation of w, A represents the matrix of moments computed from the propagator $\mathcal{D}$ and b is the vector of values of f in the mesh of evaluation points distributed within the control regions together with the prescribed two far field directions. The matrix A is not invertible in most cases, thus the linear system (24) is solved using a regularization routine to minimize the sum of squared residuals. Following the strategy in [16, 20, 21], the unknown coefficients in w_d are obtained by using Tikhonov regularization, and can be compactly written as

$${\widehat{\mathbf{w}}}_d=\underset{{\mathbf{w}}_d\in \mathrm{\partial }{D}_a^{\text{'}}}{\arg \min }\left[{\eta }_l{\Vert \mathbf{A}{\mathbf{w}}_d-\mathbf{b}\Vert }_{L^2\left(\mathrm{\partial }{W}_l\right)}^2+\alpha {\Vert {\mathbf{w}}_d\Vert }_{L^2\left(\mathrm{\partial }{D}_a^{\text{'}}\right)}^2\right] ,$$(25)

where l = 1 or 2 denoting the index of the control region. α is the regularization parameter representing the penalty weight for the power required by the solution. The optimal α is determined by the Morozov discrepancy principle [39, 40]. The unknown discrete coefficients in w_d are taken to be the Tikhonov solution,

$${\mathbf{w}}_d=(\alpha \mathbf{I}+{\mathbf{A}}^{\mathrm{*}}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\mathrm{*}}\mathbf{b},$$(26)

where I is the identity matrix and A^* is the complex conjugate transpose of A.

Remark 2.2

In (25), the weighting coefficient η_l can be applied to change the contributions of two data misfit in near region and far direction. The weighting coefficient determines the control effort made by the source in each region. For example, if η₁ → ∞, the data misfit in the near region 1 will be minimized to approach 0 such that the total misfit is less than a desired threshold. Note that the weighting coefficient also plays an important role to the complexity of the boundary input.

To estimate the power on the actual source D_a, the average radiated power P_ave and the stored energy P_stor are computed as,

$$P_{\mathrm{ave}}=\frac{1}{2}{\int }_{\mathrm{\partial }}\mathrm{}\mathrm{Re}\left[u^{\mathrm{*}}\left(\nabla u·\mathbf{n}\right)\right]\mathrm{ }\mathrm{d}S\hspace{1em}\mathrm{and},$$(27)

$$P_{\mathrm{stor}}=\frac{1}{2}{\int }_{\mathrm{\partial }}\mathrm{}\mathrm{Im}\left[u^{\mathrm{*}}\right(\nabla u·\mathbf{n}\left)\right]\mathrm{ }\mathrm{d}S,$$(28)

where u^* denotes the complex conjugate and B_R is some sphere of radius R containing the actual source D_a. In our analyses, the calculated power is expressed in dB relative to a reference level of 10⁻¹² W. Finally, the sought boundary input: either normal velocity v_n or pressure p on the actual source is obtained from (3) and (4), respectively.

3 Numerical results in free space

In this section, we present several relevant numerical simulations to support the above mentioned theoretical framework. We start from a simplified geometric configuration as shown in Figure 3a with one near control region D₁ and one far-field direction x₁ which is exactly behind the near control. Then, we extend our numerical study into a multiple-region regime with two near field control regions D₁ and D₂ and two far field directions x₁ behind D₁ and x₂ behind D₂ as sketched in Figure 3b. The source and control regions are in free space (with medium parameters c = 343 m/s and ρ = 1.225 kg/m³). Throughout this section, the fictitious source region is the sphere of radius 0.2 m centered at the origin. In general, the actual source D_a can be arbitrarily shaped as long as it is Lipschitz and compactly embeds $D_a^{\text{'}}$. In our simulations, for exemplification, we assume the actual source to be the sphere of radius 0.22 m centered at the origin. In the spirit of [35], the section starts with Sections 3.1 and 3.2 discussing the performance of our strategy, which is quantified by the relative or absolute error in control regions, in each of the above mentioned configurations and then continues with Section 3.3 where we present a detailed sensitivity analysis for free space.

Figure 3 Sketch of the top view of the problem geometry showing the near control(s) and the far field direction(s). (a) One near control and one far field direction. (b) Two near controls and two far field directions.

3.1 A null in the near control region and a non-zero far field pattern

In this subsection we show the performance of our scheme in creating a null in D₁ and a given pattern f_∞,1 in the prescribed far field direction x₁ (see Fig. 3a for one possible configuration). In our simulation, we set the wavenumber to k = 10 and consider D₁ ⋐ W₁ where W₁ is an annular sector given in the spherical coordinates (with respect to the origin) by,

$$W_1=\left\{(r,\theta ,\varphi )\mathrm{ }:\mathrm{ }r\in \left[\mathrm{0.4,0.7}\right],\hspace{1em}\theta \in \left[\frac{\pi }{4},\frac{3\pi }{4}\right],\hspace{1em}\varphi \in \left[\frac{3\pi }{4},\frac{5\pi }{4}\right]\right\}.$$(29)

The far-field direction is exactly behind the near control, i.e., ${\mathbf{x}}_1=(r,\theta ,\varphi )=(r,\frac{\pi }{2},\pi )$, for large r. The desired field in region D₁ is f₁ = 0 while the desired far field pattern in direction x₁ is given by f_∞,1 = 0.01 + i · 0.02 and $i=\sqrt{-1}$. The simulation results are shown in Figures 4 and 5.

Figure 4 Generated field on the plane z = 0.

Figure 5 Boundary input vn on the actual source ∂Da.

The generated field on a 2D cross section through z = 0 around the control region is shown in Figure 4. This shows that the field inside the control region has amplitude not exceeding −120 dB. In the linear scale, the generated field’s maximum pointwise absolute value is less than 1.6 × 10⁻⁸ while its L² norm is just about 6.69 × 10⁻⁸. As a numerical stability check, these values are computed using points slightly off from the mesh points used in the collocation scheme. To describe the accuracy of the scheme, we use the pointwise measure of error e_i defined as,

$$e_i=\{\begin{array}{c}\begin{array}{cc}\frac{|u_i-f_i|}{\left|f_i\right|}& \mathrm{if}\mathrm{ }{f}_i\ne 0,\end{array}\\ \begin{array}{cc}|u_i-f_i|& \mathrm{if}\mathrm{ }{f}_i=0,\end{array}\end{array}$$(30)

where u = Aw_d is the generated field and f_i is the prescribed value in the ith evaluation point. In the far field direction x₁, the maximum pointwise error of the real and imaginary parts are 4.79 × 10⁻⁸ and 1.68 × 10⁻⁸, respectively suggesting good control accuracy.

Remark 3.1

In Figure 4, once a scatterer besides its reflection appears in the energy radiation direction, two approaches can be applied to treat the reflections. The first method modifies the surface input of the source accounting for nulls around the scatterer, i.e., creating nulls around the scatterer such that the energy would be by design radiated away from the scatterer. Theoretically, it is possible in far field or for small scaterrers in near field. If the scatterer is large or there is a scattering boundary, the second method which incorporates the scattered field in the total propagator has to be used. This requires a different mathematical formulation and is not considered in this paper.

Figure 5 shows the normal velocity v_n on the source ∂D_a in a rectangular (θ, ϕ)-plot. The pointwise amplitudes of v_n are quite small and within order 10⁻². The average power and stored energy in the actual source are 7.1 dB and 7.8 dB, respectively, which implies that it is feasible for physical implementation, insofar as power budget is concerned.

Remark 3.2

The boundary input data shown in Figure 5 can approximate the field patterns in the exterior regions and far directions. In further numerical simulations [41], our continuous sources are proven to be able to achieve controls at regions that are much (10 wavelengths) further away from the source and that are much (15 times) larger than the source. This suggests that our approach may be useful to push the current understanding about the physical limits.

3.2 Two near control regions and two far field directions

In this subsection we show the performance of our scheme in creating null fields in D₁ and D₂ while approximating two distinct prescribed patterns f_∞,1, f_∞,2 in far-field directions x₁ and x₂, respectively (see Fig. 3b for one possible configuration). One possible application of this configuration would be in developing a strategy to establish and maintain communication in several given far-field directions while avoiding near field interference located in D₁ and D₂. In this test, k = 10, D₁ ⋐ W₁ with W₁ defined in (29) while D₂ ⋐ W₂ with W₂ given in spherical coordinates as,

$$W_2=\left\{(r,\theta ,\varphi )\mathrm{ }:r\in \left[\mathrm{1,1.2}\right],\hspace{1em}\theta \in \left[\frac{\pi }{4},\frac{3\pi }{4}\right],\hspace{1em}\varphi \in \left[\frac{-\pi }{4},\frac{\pi }{4}\right]\right\}.$$(31)

In this geometry, ${\mathbf{x}}_{\mathbf{1}}=(r,\theta ,\varphi )=(r,\frac{\pi }{2},\pi )$ and ${\mathbf{x}}_{\mathbf{2}}=(r,\theta ,\varphi )=(r,\frac{\pi }{2},0)$, r ≫ 1 while the prescribed far-field pattern values are, f_∞,1 = 0.01 + i·0.02 in direction x₁ and f_∞,2 = 0.05 + i · 0.03 in direction x₂. The simulation results are shown in Figures 6 and 7. In Figure 6, we see that the absolute errors in the two near controls are within order 10⁻⁷. In the far directions x₁ and x₂, the supremum relative errors of the real parts are 3.29 × 10⁻⁷ and 6.39 × 10⁻⁸ while for the imaginary parts, the supremum errors are 8.02 × 10⁻⁸ and 1.61 × 10⁻⁷, respectively. In Figure 7 we present the normal velocity required on the actual physical source D_a. The average radiated power is 8.4 dB and the stored energy is 9.7 dB. The power budget is slightly larger than that in Section 3.1 as the source is projecting two non-zero far-field patterns in two directions.

Figure 6 Pointwise magnitudes of the generated fields on (a) D1 and (b) D2, both approximating a null field.

Figure 7 Normal velocity vn on the actual source ∂Da.

3.3 Sensitivity analysis

The aforementioned results support the analysis of [35] and show that our strategy works for each of the two configurations depicted in Figure 3. In the following tests, we aim to study the sensitivity of our strategy with respect to variations in several physically relevant parameters, such as wavenumber k, the distance between the control region and the source, the control region size and, mutual distance between the control regions (in the case of more control regions and far field directions Fig. 3b). The feasibility of the active control scheme is also discussed by looking at the overall control accuracy, power budget and its performance against process error. The geometry in the sensitivity analysis is depicted in Figure 8. To distinguish the original region (dark), the regions with modified parameters are shown in a light color. For instance, $D_1^{\mathrm{*}2}$ denotes the near control which is shifted away from the source, in which the superscript `2’ corresponds to the experiment number in Section 3.3.2.

Figure 8 Sketch of the geometry in sensitivity analysis. D1, D2 and x1, x2 are original near controls and far field directions, respectively. They are shown in dark color. The light-colored regions with $D_1^{\mathrm{*}n}$, where n = 2, 3 and 4, are corresponding to the experiments in Sections 3.3.2 to 3.3.4.

3.3.1 Varying the wavenumber k

We start with the initial geometry in Figure 3a, i.e., only one near control and a far-field direction. The prescribed field in D₁ is zero and the far-field pattern is a given non-zero complex number. In this simulation we let the wavenumber vary from 1 to 31. In every single simulation, we keep the geometry fixed and only change the wavenumber. The simulation results are shown in Figure 9. The first plot in Figure 9 shows the supremum pointwise absolute error in D₁ and the relative error in the far field direction x₁. The behavior is less optimal in the low frequency case since the condition number of the associated matrix grows large with decrease in frequency while our numerics suggest the performance of the strategy gets worse for larger wavenumbers (k > 200). The second plot is the process error, which is the supremum pointwise error when the computed feeding current w_d is contaminated with small noise. Mathematically, the process error can be obtained by replacing u with u_δ in (30), where u_δ = Aw_dδ and δ = 0.0001 is the noise threshold. In this paper, we consider a random Gaussian noise such that w_dδ = w_d·(1 + δ·$\Vert {\mathbf{w}}_d\Vert$₂·$\mathcal{R}$) and $\mathcal{R}\sim \mathcal{N}$(0,1). The plot of the process error in Figure 9 indicates that the system is capable of overcoming feeding-current noise of a certain level. The power budget and L² norm of v_n show that the control effort decreases as the operating frequency increases.

Figure 9 Results showing the control accuracy and the power budget varying with k. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.

3.3.2 Varying the distance between the near control region and the source

The following results show the effect of variations in the distance between the near control region and the source. We still use the initial model in Figure 3a. The near control D₁ is shifted further away from the source ($D_1^{\mathrm{*}2}$ in Fig. 8a) while all other parameters are fixed. The results are depicted in Figure 10. We notice that the control accuracy in the near region keeps on improving when the near control is moved further away from the source, while the control accuracy in the far direction is converging to a certain value. This indicates that the generated field on D₁ becomes closer to the null field as it is moved further away from the source. This may be desired for scenarios when D₁ is an acoustic scatterer. The power budget on the source follows a similar trend as the control accuracy.

Figure 10 Results showing the control accuracy and the power budget varying with mutual distance between D1 and $D_a^{\text{'}}$. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. L2 norm of normal velocity vn on Da.

3.3.3 Varying the near control region size

Now we consider the behavior of the control accuracy and the power budget with respect to incremental increase in the outer radius of the near control region D₁ ($D_1^{\mathrm{*}3}$ in Fig. 8a) with all the other parameters kept fixed. The results are shown in Figure 11. Notice that the control accuracy and power budget are slightly oscillating in the entire range of near control size. The results indicate a good performance for larger obstacles.

Figure 11 Results showing the control accuracy and the power budget varying with the size of the near control D1. The size of D1 is the difference between the outer and inner radii, i.e., the thickness of the sectorial region. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.

3.3.4 Varying the mutual distance between the near control regions

In this sensitivity test, we consider two near control regions and vary the mutual distance between them. In Figure 8b, D₁ is fixed and we rotate D₂ around D_a to obtain a new secondary control region $D_2^{\mathrm{*}4}$. For a fair comparison, the far-field directions are kept exactly behind the controls, i.e., the second far field direction is also rotated (see ${\mathbf{x}}_2^{\mathrm{*}4}$). The mutual distance between the two near controls is determined by the mutual angle ϕ shown in Figure 8b. The control accuracy and power budget are recorded as the angle ϕ is varied from 3.6° to 356.4°.The results are shown in Figure 12. We find that the relative error is less than 10⁻⁵ and the power budget are within 12 dB if ϕ is in the range of [30°, 330°]. If ϕ is out of this range, i.e., the near controls are too close to each other, the relative error and source power become excessive. In the lower and upper bounds of the ϕ-range, the relative error in the far-field patterns can be as high as 15%. This indicates that the control regions cannot be too close, otherwise the accurate control effects are not guaranteed. This makes sense as it is hard to generate two very different patterns in two far field direction which are very close due to the continuity of Helmholtz fields.

Figure 12 Results showing the control accuracy and the power budget varying with the mutual distance between near controls D1 and D2 (given as a function of the angle ϕ between them). From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.

4 Numerical results in the homogeneous ocean

In this section, we extend our sensitivity analysis to the homogeneous ocean regime. We follow a similar procedure as in the previous section that deals with the problem in free space. We only show part of the numerical results due to the page limitation. The rest of the numerical examples are available in the supplementary document [41]. In the entire section, the fictitious source $D_a^{\text{'}}$ is the sphere of radius 0.2 m with center at (0, 0, −50) while, for exemplification, the physical source is chosen to be the concentric sphere of radius 0.22 m. The control region D₁ and the far-field direction x₁ are given by the control region and far-field direction described in the free space simulations but this time shifted 50 m downward (see Fig. 2). The depth of the ocean environment is |h| = 100 m, the speed of sound c is assumed to be 1515 m/s and the density ρ is 1020 kg/m³.

4.1 A null in the near control region and a non-zero far field pattern

We start from the initial geometry in Figure 3a but with k = 1 and recall, in the spirit of [35], the performance of our strategy. In D₁, we prescribe a null field, and the desired far-field pattern value is f_∞,1 = 0.01 + i·0.02. The results are shown in Figures 13 and 14. Figure 13 shows the generated fields on two cross sections of the ocean. The left plot (a) shows that the generated field on W₁ has pointwise amplitude less than −60 dB. In the linear scale, this corresponds to a maximum magnitude of the generated field within order 10⁻⁵. The right plot (b) shows the field on the plane y = 0, where the reflections from the ocean’s surface and bottom are observable.

Figure 13 Generated fields on cut planes: (a) z = −50 and (b) y = 0

Figure 14 Normal velocity vn on the actual source ∂Da.

In the exact direction x₁, the relative errors of the real part and imaginary part are 5.24 × 10⁻⁴ and 1.83 × 10⁻⁴, respectively. In Figure 9, we show the boundary input v_n on the actual source ∂D_a. The average power and stored energy on the actual source are 11.2 dB and 11.6 dB, respectively. This low power budget suggests the possibility of physically instantiating this source.

4.2 Sensitivity analysis

In this subsection we repeat the same experiments performed in the free space environment to explore the limitations or challenges of the active control in the homogeneous ocean environment. We only present here the effect of frequency and of the mutual distance between the near control and the source on the control accuracy and power budget. Some other experiment results are available online via the clickable https://drive.google.com/file/d/1E1ApCAUwBokeGuqiumu1GKzUIKR3u3N1/view?usp=sharing; Google Drive link [41].

4.2.1 Varying the wavenumber k

In this simulation, the wavenumber varies from 1 to 31 while the other parameters are fixed. The simulation results are shown in Figure 15. In contrast to the free space results shown in Figure 9, the relative error increases as the frequency increases. This is due to the complex model (propagator, boundary condition, and so on) used in the homogeneous ocean environment. In practical underwater acoustic applications, the underwater environment poses serious challenges which are much more complicated than that in the free space. Hence, the control effort required on the source to accurately manipulate the near field and directional far field pattern is more than that in free space. Furthermore, the active scheme in the underwater environment is more sensitive to the feeding noise. For the smaller noise threshold (δ = 10⁻⁸) the process error in the far direction is kept smaller than 10%. This implies that even a low noise level in the feeding network may have a major effect on the radiated pattern in a far field direction. However, the system still maintains a good performance in the near control. The overall trend of the power budget is also increasing as the wavenumber increases. It is noted that the spikes on the curves are due to the wavenumber being close to resonance frequencies.

Figure 15 Results showing the control accuracy and the power budget varying with wavenumber k. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.

4.2.2 Varying the distance between the near control region and the source

In this simulation, we test the sensitivity of our strategy as D₁ is moved away from the source as shown in Figure 8. The initial near control region is D₁ and it is incrementally pushed away from the source to obtain a new control region $D_1^{\mathrm{*}2}$. The results are shown in Figure 16. We observe that the overall performance of control accuracy and power budget is better when the near control is further away from the source, save some spikes due to resonances. The reason accounting for this trend is that the near control, where we want a field with low amplitude (as it may be interpreted as an obstacle or scatterer for some applications), is pushed away and hence the source effort is eased.

Figure 16 Results showing the control accuracy and the power budget varying with the distance between the source and the control region. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.

5 Conclusions

In this paper, the feasibility of the active manipulation of Helmholtz fields both in free space and in a homogeneous ocean of constant depth is presented. We build up on our previous works and demonstrated the possibility of of having a feasible characterization of an active source (modeled as surface pressure or surface normal velocity) such that it is capable of approximating a priori given field some near control regions while simultaneously projecting desired patterns in several far-field directions. By a feasible source, we mean a source with low power requirement, for instance 30 dB, which can radiate the prescribed fields within some accuracy thresholds. We showed the control accuracy and the power budget of the proposed active control mechanism for each environment. Then we explored the behavior of physically relevant parameters (power budget and control accuracy) with respect to variations in the frequency, outward shift, the outer radius of the near control and the mutual distance between near controls.

In our simulations, we considered the two initial models shown in Figure 3. The first one contains one near control and one far-field direction and the second one has two near control regions and two far-field directions. The far-field directions are placed exactly behind the near control regions in each of the two models. In this paper, we only show the cases in which the near control regions are prescribed to be have a null field and the given far-field directions to have a non-zero pattern.

In free space, the control accuracy is within order 10⁻⁸ both in the near control and in the far field direction. In the geometry shown in Figure 3a, the operating frequency is first varied while all other parameters are fixed. The frequency is swept from 54.59 Hz to 1.69 kHz (k is from 1 to 31). Moreover, the power requirement is kept at low levels which suggests the feasibility of a physical implementation of the calculated source. In the second simulation, the near control region is moved outward. We noticed that the control accuracy in the near region is continuously improving as the near control region is pushed further away from the active source. However, accuracy converges to a certain level in the far-field direction. Next, we vary the outer radius of the near control region to explore the effect of the near region’s size on the control accuracy and power budget. The simulation results show that the active control scheme works in the entire range of values of the outer radii. We also considered the geometry in Figure 3b. Here, we rotate the second near control region together with the far-field direction behind it. We find that the control accuracy and power budget don’t change significantly if the two near control regions are separated by an angle ϕ ∈ [30°, 300°]. However, outside this range the control performance is gradually degrading.

We then extended our sensitivity analysis into the case of a homogeneous ocean of constant depth. The Green’s function, corresponding to a pressure-release surface and a totally reflecting bottom, is expressed using the normal mode representation. Accordingly, the far-field pattern propagator is defined. In this case, the control accuracy is within orders 10⁻⁵ and 10⁻⁴ in the near control and in the far direction, respectively. In the first sensitivity test we varied the wavenumber from 1 to 31. The results show that the overall performance of the scheme, with respect to the control accuracy and the power budget decreases with increase in frequency. These suggest that the control scheme is more suitable to low frequencies in the homogeneous ocean environment. In the second test, we move the near control further away from the source. The results are similar to that in the free space. More sensitivity test results are available in the supplementary material [41].

Conflict of interest

Author declared no conflict of interests.

Acknowledgments

This research has been supported in part by the Army Research Office under award number W911NF-17-1-0478 and in part by the National Science Foundation under award 1801925.

References

All Figures

	Figure 1 Sketch of the problem geometry showing the near controls D1, D2 and the far field directions x1 and x2 in free space.
In the text

	Figure 2 Sketch of the problem geometry showing the near controls D1, D2 and the far field directions x1 and x2 in homogeneous ocean with constant depth.
In the text

	Figure 3 Sketch of the top view of the problem geometry showing the near control(s) and the far field direction(s). (a) One near control and one far field direction. (b) Two near controls and two far field directions.
In the text

	Figure 4 Generated field on the plane z = 0.
In the text

	Figure 5 Boundary input vn on the actual source ∂Da.
In the text

	Figure 6 Pointwise magnitudes of the generated fields on (a) D1 and (b) D2, both approximating a null field.
In the text

	Figure 7 Normal velocity vn on the actual source ∂Da.
In the text

	Figure 8 Sketch of the geometry in sensitivity analysis. D1, D2 and x1, x2 are original near controls and far field directions, respectively. They are shown in dark color. The light-colored regions with $D_1^{\mathrm{*}n}$, where n = 2, 3 and 4, are corresponding to the experiments in Sections 3.3.2 to 3.3.4.
In the text

	Figure 9 Results showing the control accuracy and the power budget varying with k. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.
In the text

	Figure 10 Results showing the control accuracy and the power budget varying with mutual distance between D1 and $D_a^{\text{'}}$. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. L2 norm of normal velocity vn on Da.
In the text

	Figure 11 Results showing the control accuracy and the power budget varying with the size of the near control D1. The size of D1 is the difference between the outer and inner radii, i.e., the thickness of the sectorial region. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.
In the text

	Figure 12 Results showing the control accuracy and the power budget varying with the mutual distance between near controls D1 and D2 (given as a function of the angle ϕ between them). From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.
In the text

	Figure 13 Generated fields on cut planes: (a) z = −50 and (b) y = 0
In the text

	Figure 14 Normal velocity vn on the actual source ∂Da.
In the text

	Figure 15 Results showing the control accuracy and the power budget varying with wavenumber k. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.
In the text

	Figure 16 Results showing the control accuracy and the power budget varying with the distance between the source and the control region. From left to right, (1) Supremum error. (2) Process error. (3) Power budget on Da. (4) L2 norm of normal velocity vn on Da.
In the text