# Isogeometric simulation of acoustic radiation

In this paper we discuss the numerical solution on a simple 2D domain of the Helmoltz equation with mixed boundary conditions. The so called radiation problem depends on the wavenumber constant parameter k and it is inspired here by medical applications, where a transducer emits a pulse at a given frequency. This problem has been successfully solved in the past with the classical Finite Element Method (FEM) for relative small values of k. But in modern applications the values of k can be of order of thousands and FEM faces up several numerical difficulties. To overcome these difficulties we solve the radiation problem using the Isogeometric Analysis (IgA), a kind of generalization of FEM. Starting with the variational formulation of the radiation problem, we show with details how to apply the isogeometric approach in order to compute the coefficients of the approximated solution of radiation problem in terms of the B-spline basis functions. Our implementation of IgA using GeoPDEs software shows that isogeometric approach is superior than FEM, since it is able to reduce substantially the pollution error, especially for high values of k, producing additionally smoother solutions which depend on less degrees of freedom.

## Authors

• 4 publications
• 2 publications
• 3 publications
• 2 publications
• 2 publications
01/21/2020

### Isogeometric solution of Helmholtz equation with Dirichlet boundary condition: numerical experiences

In this paper we use the Isogeometric method to solve the Helmholtz equa...
09/13/2019

### Finite element modeling of micropolar-based phononic crystals

The performance of a Cosserat/micropolar solid as a numerical vehicle to...
01/24/2022

### The Helmholtz boundary element method does not suffer from the pollution effect

In d dimensions, approximating an arbitrary function oscillating with fr...
01/11/2021

### Adaptive modelling of variably saturated seepage problems

07/19/2018

In this paper, first we present the variational formulation for a second...
02/21/2022

### Isogeometric Analysis of Bound States of a Quantum Three-Body Problem in 1D

In this paper, we initiate the study of isogeometric analysis (IGA) of a...
02/14/2022

### A simple proof that the hp-FEM does not suffer from the pollution effect for the constant-coefficient full-space Helmholtz equation

In d dimensions, approximating an arbitrary function oscillating with fr...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Wave problems have been intensively studied since they are relevant in multiple fields. The solution of wave equation is usually written as the product of a function of time and a function which only depends on spacial variables. In acoustic problems, for instance, the time function is chosen as , where is the angular frequency of the propagating wave and i is the imaginary unit. With this time harmonic dependence, the wave equation is reduced to the Helmholtz equation

 −△u(x,y)−k2u(x,y)=0

where is the number of waves per unit of distance, called wavenumber, and is the wavelength.

Helmholtz equation is very important in acoustic applications, including seismic wave propagation, acoustic noise control, non destructive testing and medical ultrasound. In particular, therapeutic applications of ultrasound involve focused beams directing the ultrasound energy into the tissue region that needs the treatment. Currently, High Intensity Focused Ultrasound (HIFU) therapy method is known as one of the most advanced surgical [23] and also physio-therapeutical techniques [18]. In most clinical applications, HIFU transducers are excited at a single frequency in the range MHz. From the mathematical point of view, Pennes’ bioheat equation [30] is used to model thermal diffusion effects of HIFU. This equation relates the temperature distribution in time and space with the absorbed ultrasound energy, which is computed from the acoustic pressure field solution of the Helmholtz equation.

The numerical solution of Helmholtz equation is in general a challenge. When the wavenumber is small, it can be handled using low order Finite Element Method (FEM). But the design of robust and efficient numerical algorithms for high values of is difficult. In practice, many numerical difficulties appear. First, since the function oscillates on a scale of , to obtain an accurate approximation of using finite elements of degree and increasing , it is necessary to request that the total number of degrees of freedom to be proportional to [17], or to choose a mesh size such that is constant and sufficiently small [28]. It means that very large linear systems have to be solved, with high computational cost. Moreover, it is known [1], [2], [22], that even if we use a big number of degrees of freedom, the errors of continuous Galerkin finite element approximations increases when becomes larger. In the literature, this non-robust behavior with respect to , is known as the pollution effect. To reduce the pollution several authors have proposed to enrich the basis of the finite element space with wave-like functions depending on the wavenumber. One of the first steps in this direction is the Partition of Unity Finite Element Method [27], [29].

The standard variational formulation of the Helmholtz equation is sign-indefinite (i.e.not coercive). Hence, another difficulty for the numerical solution of the Helmholtz equation is that for sufficiently large, the coefficient matrix is indefinite and non-normal. As a consequence, iterative methods to solve the corresponding linear systems behave extremely bad if the system is not preconditioned [14], [11]. To face this problem researches have proposed several preconditioners, such as multigrid methods with Krylov smoothers, domain decomposition, and complex shifted Laplacian preconditioner. The last one was introduced in [13] and further developed and successfully generalized in [32], [17] and [33].

Dealing with wave problems, the small discrepancies between the boundary of the mesh constructed by FEM and the boundary physical domain , can significantly increase the error of the FEM approximated solution [29]. This is more evident in 3D industrial applications, where the surface of the physical domain is usually represented in terms of Nonuniform Rational B-spline functions (NURBS) [31]

. Since B-spline spaces include as a particular case the piecewise polynomial spaces commonly used in FEM, it was natural to think of the possibility of writing the approximated solution of the partial differential equation (PDE) in terms of the B-spline basis functions. This idea led to the emergence of the Isogeometric Analysis (IgA), introduced by Hughes et al. in 2005

[20], as a modern method to solve PDE. IgA uses B-spline functions to parametrize the geometry and as shape functions to approximate the solution of the PDE. In this sense, it combines the variational techniques of isoparametric FEM, with the classical functions in computer design systems. IgA and FEM are based on the same principle, the Galerkin method, but IgA approach has a very important advantage: B-spline basis functions may be constructed to have high smoothness. This is crucial in problems with smooth solutions, where improved accuracy per degree of freedom is obtained in comparison with the classical FEM. It explains the wide range of applications solved successfully in the last years with IgA approach, see for instance [7], [4], [34], [5], [19].

### 1.1 Related work

The literature dealing with different aspects of the solution of Helmholtz equation with IgA is very recent [24],[8],[10], [11] and [12].

In [24] the performance of IgA to solve exterior scattering problems is investigated, using an absorbing boundary condition on a fictitious boundary to truncate the infinite space. It is shown that IgA is a robust approach to reduce the effects of the pollution error and therefore it is a promising tool to solve high-frequency acoustic problems. In [8] IgA is used to solve Helmholtz equation with several boundary conditions in 2D regions. The results of a convergence study are presented confirming that IgA outperforms FEM for similar degrees of freedom, specially when the frequency of the waves increases. In [10] the effect of higher continuity of B-spline basis function on the pollution error is studied. The conclusion is that the pollution is improved with IgA compared to classical FEM. Moreover, it is shown that partition of unity isogeometric analysis (PUIgA) is suitable for wave problems, since enrichment eliminates the need of domain re-meshing at higher frequencies.

In [11] the Helmholtz equation with Robin boundary condition is tackled using IgA. GMRES method for solving the linear system resulting from IgA is investigated, including the use of preconditioners such as ILU with a complex shift and complex shifted Laplace. The study in [11] concludes that, for all wavenumbers, GMRES converges at a fewer iterations with IgA compared to FEM. Moreover, the pollution error is significantly reduced with IgA, even when it is not completely eliminated. In a very recent paper [12], the focus of the research is on the numerical solution of the linear system derived from IgA discretization of Helmholtz equation. The system is solved with GMRES and its convergency is accelerated using a deflation technique, combined with the approximated computation of the inverse of the CSLP with a geometric multigrid method. Numerical results for one a two dimensional problems are shown, confirming scalable convergence with respect to the wavenumber and the order of the B-spline basis functions.

### 1.2 Our contribution

The main contribution of this paper is the application of IgA to the solution of a radiation problem, mathematically modeled with a 2D Helmholtz equation with mixed boundary conditions. All the details concerning the application of IgA method are presented, including the computation of the matrix and the right hand side vector of the linear system, whose solution provides the B-spline coefficients of the approximated solution. With the open source software GeoPDEs

[15] we have implemented our in-house code to solve the radiation problem using IgA. The results of the numerical implementation with IgA are compared with the solution using classic FEM, for relative high values of the wavenumber. We show the superiority of the performance of IgA approach by means of several experiments, which confirm that using less degrees of freedom smoother approximated solutions are obtained with a substantially reduced pollution.

## 2 Physical problem and variational formulation

In this paper we are interested in acoustic wave applications. Under the assumption that the acoustic wave propagation is linear and also that the amplitude of shear waves in the media are much smaller than the amplitude of the pressure waves, nonlinear effects and shear waves may be neglected. In consequence, the acoustic wave pressure is a complex function solution of the Helmholtz equation and the wavenumber is a positive and real number.

Inspired by the experiments to measure focused ultrasound induced heating in a tissue phantom [3], [26] we consider a 2D axial symmetry geometry : the semicircle of radius and center on the origin of coordinates. Moreover, we assume that a transducer of aperture , with , and flat geometry is located at the origin, see Figure 1. The transducer emits a piston-like pulse of frequency , with speed and constant amplitude equal to . Hence, Dirichlet boundary condition is imposed on . Additionally, boundary is simulated as an acoustically rigid wall by setting the normal velocity equal to zero. Dirichlet and Neumann boundary conditions are known in the literature as rigid and free baffle respectively. Finally, assuming that we require that the boundary has the same acoustic impedance as the media to avoid wave reflections. This Robin condition is referred to as impedance boundary condition.

In the rest of the paper we call radiation problem to the solution of equation

 −△u(x,y)−k2u(x,y)=0,(x,y)∈Ω (1)

with mixed boundary conditions

 u(x,y) = ConΓD (2) ∂u(x,y)∂→n = 0onΓN (3) ∂u(x,y)∂→n+iku(x,y) = 0onΓR (4)

where denotes the normal vector to the boundary or .

Near the transducer, in the near field area, there are significant fluctuations in the ultrasound intensity. However, from some point on the pressure waves form a relatively uniform front that spreads out in a pattern originating from the center of the transducer. This area is called the far field and it is important in applications, since optimal detection occur at the start of far field, where the sound wave is well behaved and attains its maximum strength. The near field length defines the transition point between the near field and the far field. This point, sometimes referred to as the “natural focus”, can be calculated as

 Nf=a2λ (5)

where is the wavelength.

### 2.2 Variational formulation

In this section we obtain variational formulation of the radiation problem. Denote by the Hilbert space

 V0={v:Ω→C,v∈H1(Ω),v(x,y)=0% for(x,y)∈ΓD} (6)

with the norm

 ∥v∥2V0:=∥∇v∥2L2(Ω)=∫∫Ω∇vt∇¯¯¯vdΩ (7)

where is the complex conjugate of and is the gradient vector of the complex function . To obtain the variational formulation we multiply (1) by with and integrate on . Applying Green’s first identity and imposing the mixed boundary conditions we arrive to the variational formulation,

 Findu∈H1(Ω)withu=ConΓD, such thata(u,v)=0for allv∈V0 (8)

where the sesquilinear form is given by,

 a(u,v) = ∫∫Ω(∇u(x,y)t∇¯¯¯v(x,y)−k2u(x,y)¯¯¯v(x,y))dΩ+ik∫ΓRu(x,y)¯¯¯v(x,y)ds (9)

The solution of the variational problem (8) is called weak solution. We prove that it exits considering two auxiliary problems and writing the weak solution as the sum of the solutions of the auxiliary problems.

###### Lemma 1

For each there exits a weak solution of the auxiliary problem

 −△u0(x,y)−k2u0(x,y)=f(x,y),(x,y)∈Ω (10)

with homogeneous mixed boundary conditions

 (11)

Proof (main steps): A weak solution of the auxiliary problem is a solution of the variational problem

 Findu0∈V0, such thata(u0,v)=⟨f,v⟩for allv∈V0 (12)

where is the scalar product in and is given by (9). The form in (9) is sesquilinear and continuous. Continuity follows from the Cauchy-Schwarz inequality and the continuity of the trace map , since involves an integral over with . Although is not coercive, it satisfies a Gårding inequality . Indeed,

 Re(a(v,v))+k2∥v∥2L2(Ω)=∥∇v∥2L2(Ω)=∥v∥2V0.

Hence Fredholm theory can be applied, to show that a solution to the variational problem (12) exists since the homogeneous adjoint problem has an unique solution. Details of the proof can be found in [21].

The next result shows that a solution of (8) may be easily constructed in terms of a new function and the weak solution of the problem (10) with .

###### Proposition 1

Let be a function satisfying the conditions

 uC(x,y)=ConΓD,uC(x,y)=0onΓR,∇uC(x,y)=0onΩ (13)

Denote by the solution of (12) with . Then, the function

 u=u0+uC (14)

is a weak solution of (1) with boundary conditions (2)-(4), i.e a solution of (8).

Proof
For we have . Furthermore, for all , from (12) we get, . Moreover, from (9) and (13) it holds, . Finally, since is a sesquilinear form, from the previous equalities and (14) we obtain, for all , i.e. is a solution of (8).

Observation: function satisfying conditions (13) is not unique. In consequence, the weak solution of the radiation problem is not unique either.

## 3 Galerkin method with isogeometric approach

The Galerkin method solves the variational problem assuming that the approximated solution belongs to a finite-dimensional subspace of . In the classical FEM, consists of piecewise polynomials functions with global continuity. In the isogeometric approach [9],

is generated by tensor product NURBS functions with higher global continuity. Moreover, the physical domain

is previously parametrized by a smooth function

 F(ξ,η):^Ω⟶Ω

defined on the unit square and with piecewise smooth inverse.

### 3.1 Parametrization of the domain

In this paper the map transforming in is defined subdividing the semicircle in 3 circular sectors and with , see Figure 2. We assume that “left” and “right” curves and respectively, have both the same arc length and define an angle with the axis, with . These curves can be written exactly as quadratic rational B-spline curves for the same sequence of knots . Similarly, the “top” curve can be expressed as a quadratic rational B-spline curve for a sequence of knots . Denote by the segment of line passing through the points and . Elevating the degree of , it can be represented also as a quadratic rational B-spline curve for the sequence of knots .

The map is computed as the bilinearly blended Coon’s patch [31]interpolating the curves ,, and . It can be written as

 F(ξ,η)=nF∑i=1mF∑j=1wi,jPi,jR3i,τξ(ξ)R3j,τη(η) (15)

where and are respectively the -th and -th rational quadratic B-splines of order , for the knots and , and are the weights.

The map satisfies the interpolation conditions

 F(ξ,0) = cb(ξ),F(ξ,1)=ct(ξ),0≤ξ≤1 F(0,η) = cl(η),F(1,η)=cr(η),0≤η≤1

The points are the vertices of the control mesh of the map . For , the curves can be represented as rational quadratic B-spline curves with knots , in consequence . Similarly, the curve can be written as a rational quadratic B-spline curve with knots , hence . The 9 control points and the control mesh of are shown in second row of Figure 7 (left).

The quality of the parametrization of the physical domain influences the precision of the solution computed with isogeometric approach [25]. In this sense, a good uniformity and orthogonality of the isoparametric curves of is desirable. In the literature the quality of the parametrization in the point is measured computing the mean ratio Jacobian given by,

 Jr(ξ,η)=2detJF(ξ,η)∥Fξ(ξ,η)∥22+∥Fη(ξ,η)∥22 (16)

where and are the tangent vectors to the isoparametric curves, denotes the Euclidean norm and denotes the Jacobian matrix of the parametrization,

 JF(ξ,η)=(xξxηyξyη)

If the map is injective, then does not changes of sign. Assuming that it is positive it holds that . A value of equal to 1 at a point indicates that the isoparametric curves are orthogonal at and the map produces the same length distortion at in both parametric directions and . In the center of Figure 7 we show the mesh in with vertices computed as the image by of the vertices of a rectangular mesh in . For the corresponding parametrization preserves the geometry of the quadrilateral in is almost everywhere, except in the areas near the two points subdividing the semicircle in three curves. This is observed in the center of the second row of Figure 7. On the right of this figure we show a color map, where colors correspond to the values of according to (16). A yellow color indicates that . Hence, the distortion introduced by the parametrization for is small almost everywhere, except on the two blue areas (values of close to ), where the distortion produced by the parametrization is higher. Observe that blue areas are not contained in the region , where the highest pressure values are located, see Figure 3 (left).

### 3.2 Variational formulation on [0,1]×[0,1]

With the help of the parametrization the double integral in (9) may be transformed into an integral over according to the integration rule

 ∫∫Ωh(x,y)dΩ=∫10∫10h(F(ξ,η))|detJF(ξ,η)|dξdη

Moreover, from the chain rule applied to

we know that . Hence, the double integral in (9) can be written as

 a(u,v) = ∫10∫10(JF(ξ,η)−t∇u(F(ξ,η)))t(JF(ξ,η)−t∇¯¯¯v(F(ξ,η)))|detJF(ξ,η)|dξdη (17) − k2∫10∫10u(F(ξ,η))¯¯¯v(F(ξ,η))|detJF(ξ,η)|dξdη+ik∫ΓRu(x,y)¯¯¯v(x,y)ds

To obtain a formulation in the parametric domain it is necessary to rewrite the last integral in (17). Recall that the parametrization has been constructed in such a way that and has been subdivided in three consecutive boundaries with

 ΓR1 = F(0,η)=cl(η),0≤η≤1 (18) ΓR2 = F(ξ,1)=ct(ξ),0≤ξ≤1 (19) ΓR3 = F(1,η)=cr(η),0≤η≤1 (20)

Taking into account that , from (18),(19) and (20) we obtain an expression for that only depends on variables and . Substituting this expression in (17) we obtain finally that can be written in the parametric domain as,

 a(u,v) = ∫10∫10∇u(F(ξ,η))t(JF(ξ,η)tJF(ξ,η))−1∇¯¯¯v(F(ξ,η))|detJF(ξ,η)|dξdη − k2∫10∫10u(F(ξ,η))¯¯¯v(F(ξ,η))|detJF(ξ,η)|dξdη + ik∫10u(F(0,η))¯¯¯v(F(0,η))⎛⎝(∂x(0,η)∂η)2+(∂y(0,η)∂η)2⎞⎠1/2dη + + ik∫10u(F(1,η))¯¯¯v(F(1,η))⎛⎝(∂x(1,η)∂η)2+(∂y(1,η)∂η)2⎞⎠1/2dη

### 3.3 Galerkin method with B-spline functions

In the isogeometric approach the Galerkin method replaces the infinite dimensional space by the space generated by tensor product B-spline functions [31]. These functions are computed using univariate B-splines in both parametric directions and . To define the B-splines of order (degree ) in the direction we need a sequence of knots . Similarly, the definition of the univariate B-splines of order (degree ) in the direction requires a sequence of knots . Given and , the knots sequences are defined in this work by

 tξ = (k1−10,...,0,ξ1,ξ2,...,ξn−k1+2,k1−11,...,1) (22) tη = (k2−10,...,0,η1,η2,...,ηm−k2+2,k2−11,...,1) (23)

where are the breakpoints in the direction and are the breakpoints in the direction . The sequences of breakpoints define a rectangular mesh in with vertices .

The B-splines functions , are a basis of the space of splines of order with knots . Similarly, the B-spline functions are a basis of the space of spline functions of order with knots . The functions,

 Bk1,k2i,j(ξ,η):=Bk1i,tξ(ξ)Bk2j,tη(η),i=1,...n,j=1,...,m (24)

are a basis of the tensor product space of splines functions of order in the direction and order in the direction . To simplify the notation, in the rest of the paper we don’t write the subindex or of the B-spline functions when it is clear from the context. Due to the assumptions on the parameterization , the functions

 ϕi,j(x,y)=(Bk1,k2i,j∘F−1)(x,y),i=1,...n,j=1,...,m (25)

are linearly independent in . Denote by the space,

 Vh=span{ϕi,j(x,y),i=1,...,n,j=1,...,m} (26)

Galerkin method computes the approximated solution of the variational problem (8) as a function in . It means that is written as,

 uh(x,y)=n∑i=1m∑j=1γi,jϕi,j(x,y) (27)

for certain unknown coefficients . In order to obtain a system of linear equations for the unknowns it is convenient to vectorize the basis functions and the corresponding coefficients in (27) introducing the change of indices

 q=q(i,j):=n(j−1)+i,i=1,...,n,j=1,...,m (28)

With this transformation we rewrite the expression (27) as

 uh(x,y)=N∑q=1αqψq(x,y) (29)

where and

 α: = (α1,...,αN)=(γ1,1,...,γn,1,γ1,2,...,γn,2,γ1,m,...,γn,m) (30) ψ(x,y): = (ψ1(x,y),...,ψN(x,y)) (31) = (ϕ1,1(x,y),...,ϕn,1(x,y),ϕ1,2(x,y),...,ϕn,2(x,y),ϕ1,m(x,y),...,ϕn,m(x,y))

Thus, in the new notation

 Vh=span{ψq(x,y),q=1,...,N}

Denote by the subspace of constituted by functions in vanishing on ,

 Vh0=span{ϕi,j(x,y), such thatϕi,j(x,y)=0for all(x,y)∈ΓD} (32)

To determine which are the functions in we have to compute the preimage by of the points and delimiting . Since with is , it is clear that there exist and both in such that,

 F(ξa−,0) = (−a,0) (33) F(ξa+,0) = (a,0) (34)

From the expression (15) it is easy to check that the solution of equations (33) and (34) are: and , both in . Assume that and satisfy

 tξi1≤ξa−

with . Since with and , for , we conclude that the B-splines not identically null on are

 ϕi,1(x,y)=Bk1,k2i,1(F−1(x,y))=Bk1i(ξ)Bk21(η),i1−k1+1≤i≤i2

Hence,

 Vh0=span{ϕi,j(x,y),1≤i≤n,2≤j≤m,% andϕi,1(x,y),i∉[i1−k1+1,i2]} (35)

Denote by the set containing the global indices, computed using (28), of functions on . Then for and . In other words, and , where is the size of . Similarly, denote by the set of global indices of functions with and let be the size of . With this notation and can be written as

 uh(x,y)=uh0(x,y)+uhg(x,y)

where and . Observe that function . Hence to obtain the Galerkin formulation we substitute in (LABEL:axieta01R) by given by (29) and by a basis function of . The result is,

 N∑q=1αq(sp,q−k2mp,q+ikep,q)=0p∈I0 (36)

where ( omitting the dependence of functions of when it is clear)

 sp,q = ∫10∫10(∇ψq)t(JFtJF)−1∇ψp|detJF|dξdη mp,q = ∫10∫10ψqψp|detJF|dξdη ep,q = ∫10ψq(F(0,η))ψp(F(0,η))⎛⎝(∂x(0,η)∂η)2+(∂y(0,η)∂η)2⎞⎠1/2dη + ∫10ψq(F(ξ,1))ψp(F(ξ,1))⎛⎝(∂x(ξ,1)∂ξ)2+(∂y(ξ,1)∂ξ)2⎞⎠1/2dξ + ∫10ψq(F(1,η))ψp(F(1,η))⎛⎝(∂x(1,η)∂η)2+(∂y(1,η)∂η)2⎞⎠1/2dη

Finally we have to impose the Dirichlet boundary condition on . For that we set,

 αq=C,q∈Ig (37)

Observe that with this assignment, for any point in it holds

 uh(˜x,0)=uh0(˜x,0)+uhg(˜x,0)=0+∑q∈Igαqψq(˜x,0)=C∑q∈Igψq(˜x,0)=C

where the last equality holds since B-spline functions satisfy the unit partition property.

Without loss of generality, assume that unknowns have been reorganized in such a way that the first unknowns correspond to indexes in and the last unknowns correspond to indexes in . Then, taking into account (37) the linear equations (36) can be written as,

 n0∑q=1αqap,q=CN∑q=n0+1ap,qp=1,...,n0 (38)

where . In the literature matrices and for are known as stiffness matrix and mass matrix respectively. Let , then system (38) can be written in matrix form as,

 A˜α=b (39)

where for and is given by

 A=S−k2M+ikE (40)

with .

This is a good place to recall that the numerical solution of the linear system (39) is a challenge. Some observations in this sense are the following. Even when matrix is symmetric, it is also indefinite and not Hermitian. Moreover, is sparse but it gets denser with the increase of the order of B-splines, since the ratio of the nonzero elements of to total number of elements of is bounded above by , with . For large values of ,