Numerical solution to the Complex 2D Helmholtz Equation based on Finite Volume Method with Impedance Boundary Conditions

Angela Handlovičová; Izabela Riečanová

doi:10.1515/phys-2016-0044

Article Open Access

Numerical solution to the Complex 2D Helmholtz Equation based on Finite Volume Method with Impedance Boundary Conditions

Angela Handlovičová and Izabela Riečanová

Published/Copyright: December 12, 2016

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 14 Issue 1

Abstract

In this paper, the numerical solution to the Helmholtz equation with impedance boundary condition, based on the Finite volume method, is discussed. We used the Robin boundary condition using exterior points. Properties of the weak solution to the Helmholtz equation and numerical solution are presented. Further the numerical experiments, comparing the numerical solution with the exact one, and the computation of the experimental order of convergence are presented.

Keywords: Helmholtz equation; Finite volume method; numerical solution; stability

PACS: 02.30.Jr; 47.11.Df

1 Introduction

Numerical methods in acoustics solve the wave equation

(1)∂2P∂t2=c2ΔP,

where P is the pressure and c is the speed of sound. This equation describes the behavior of sound, light, or water waves. In the case of time harmonic acoustic propagation and scattering [1], the pressure function is given by

(2)Px,t=ReA(x)e−iωt.

Here ω = 2πf is the angular frequency measured in rad/s, ƒ the frequency measured in Hz and Re denotes the real part. Function A is in general complex valued and is the complex acoustic pressure.

To solve (1), the method of separating the variables can be applied. Thus the time-independent form of the original equation is obtained, which is called the Helmholtz equation

(3)ΔA+k2A=0.

Here |A(x)| is the amplitude of the time harmonic pressure fluctuation at x and k is the acoustic wavenumber (number of radians per unit distance). Wavenumber is also given by the formula

(4)k=ωc.

The Helmholtz equation is related to the problems of steady-state oscillations. The unknown function A(x) is defined on a two or three dimensional domain D, where its boundary is denoted by ∂D. In our paper we focus on the most commonly relevant boundary condition [1], called impedance boundary condition of the type

(5)∂A∂n+ikβA=g.

In this boundary condition n denotes the outward normal to the boundary ∂D, and ∂∂n denotes the normal derivative. The function g on the right side of the equation is identically zero in acoustic scattering problems [1], and is nonzero for the radiation problems. It is be generally considered as the sound source. β is the relative surface admittance. In general it is a function of position on the boundary and frequency. The simplest case is when the boundary is acoustically rigid or sound hard. In this case there is no flow across ∂D and β = 0, so we obtain Neumann boundary condition. If the value is set to β = 1, it means the boundary has maximum sound absorption, so the free space is considered. The results of simulations with changing values of β and g can be seen in [2].

Our main goal is to present a numerical scheme for solving the problem (3), (5) in 2D.

2 Finite Volume Method

There are several numerical techniques for solving the Helmholtz equation. Among them we can mention the Finite element method e.g. in [3], or the Boundary element method e.g. in [4] and [5]. In this article we study numerical solution based on the Finite volume method which is an extension of the previous work [4].

We present the numerical scheme based on the Finite volume method [6]. The discretization of the domain D is the union of so called finite volumes (in 2D usually rectangles or triangles). This discretization is denoted as T_h, where the index h is connected with the size of finite volumes. In our case the domain will be a two dimensional rectangle and our finite volumes will be squares of size h. In each finite volume p ∊ T_h we have a representative point X_p in which the approximated function can be evaluated. That is why our numerical solution is a piecewise constant function, which is constant on each finite volume, and is calculated at the representative point. This point is usually chosen in the barycentre of the finite volume. Moreover we denote by E the set of all edges of each finite volume p ∊ T_h. If we have our discretization as described above, our mesh fulfils an important property

(6)Xp−Xq=dpqnpq

for both neighbouring representative points X_p and X_q. Here n_pq is the outward normal of the finite volume p to the common side with finite volume q; this side is denoted by σ_pq and d_pq is the distance between representative points X_p and X_q. An important feature of the method is the local conservativity of numerical fluxes, which means that the flux is conserved from one discretization finite volume to its neighbour.

After the discretization of the domain we have n finite volumes along one side of the rectangle domain and m finite volumes along the other, a mesh of m × n finite volumes is obtained. The particular finite volume is labelled as p and its boundary as ∂p. In the finite volume we denote the constant value of the approximated solution as u_p, and the solution in the neighbouring volume as u_q. The size of the finite volume p is denoted by m(p) and the edge σ_pq has the size denoted by m(ς_pq). We denote by N(p) the set of all neighbours of the finite volume p, which means the finite volumes that have common side with volume p.

The finite volume numerical scheme can be obtained by integrating the differential equation (3) on each finite volume. Using Green’s theorem we obtain

(7)∫pk2Adx+∫∂p∇A⋅nds=0.

We apply the approximation function as discussed above which can be denoted by

(8)uh(x)=up,x∈p.

We use this property in (7) together with the approximation of the normal derivative by a standard finite difference. This way we obtain

(9)k2upm(p)+∑q∈N(p)uq−updpqm(σpq)=0.

The solution to the Helmholtz equation is based on complex values as

(10)A=Ar+iAi.

For the approximate solution it is the same

(11)uh=uhr+iuhi,

and analogously

(12)up=upr+iupi,

where i is imaginary unit. Thus the proposed equation is valid for both real part and imaginary part. We will denote it in a similar way but with the upper indices “r” or “i”. We have

(13)k2uprm(p)+∑q∈N(p)uqr−uprdpqm(σpq)=0,k2upim(p)+∑q∈N(p)uqi−upidpqm(σpq)=0.

These equations are valid for the interior finite volumes.

We now denote by E = E_int ∪ E_ext, where E_int is the set of all interior edges and E_ext is the set of all edges of all finite volumes that belong to ∂D. Further T_h = T_h,int ∪ T_h,ext, where T_h,int is the set of all finite volumes which have all edges in E_int and T_h,ext is the set of all finite volumes that have at least one edge in E_ext. Finally by N(p)_int we denote the set of those neighbours with common side σ_pq ∊ E_int, and N(p)_ext is the set of neighbours with common side σ_pq ∊ E_ext.

From prescribed boundary condition (5), we obtained the conditions for the real and imaginary parts of the solution

(14)∂Ar∂n−kβAi=gr,∂Ai∂n+kβAr=gi

For the finite volume p ∊ T_h,ext we use the boundary conditions to approximate the numerical fluxes along the edges from T_h,ext. For this purpose we use exterior finite volumes denoted by qext, and the value of the numerical solution at these finite volumes is denoted by uqextr and uqexti, respectively. Thus for p ∊ T_h,ext we have

(15)k2uprm(p)+∑q∈N(p)intuqr−uprdpqm(σpq)+∑qext∈N(p)extuqextr−uprdpqextm(σpqext)=0,k2upim(p)+∑q∈N(p)intuqi−upidpqm(σpq)+∑qext∈N(p)extuqexti−upidpqextm(σpqext)=0.

From (14) we obtain

(16)uqextr−uprdpqext−kβuqexti+upi2=gpqr,uqexti−upidpqext+kβuqextr+upr2=gpqi,

where gpqr and gpqi are values of real and imaginary part of prescribed boundary function g evaluated on the exterior edges belonging to N(p) in the following way

(17)gpqr=1m(σpqext)∫σpqextgr(s)ds,gpqi=1m(σpqext)∫σpqextgi(s)ds.

Now, from equations (16), we eliminate the unknown values ueqxtr and ueqxti that we substitute in the initial equations (15) in the numerical fluxes along the exterior edges. For more detailed description see section 4.

This way we create a system of linear algebraic equations, in which the matrix was of order 2nm . After solving the system, both the real and imaginary part are obtained for each finite volume.

3 Properties of the Weak and Numerical Solution

Let the data in (5) fulfil the following assumptions

g = g^r + igⁱ and g^r ∈ L₂(∂D), gⁱ ∈ L₂(∂D)
β is real number

Definition

A complex valued function A = A^r + iAⁱ is a weak solution of (3)-(5) if

A^r ∊ H¹(D)) and Aⁱ ∊ H¹(D)
the following holds

(18)∫D∇A∇vdx=∫∂Dg+iβkAvds+k2∫DAvdx∀v=vr+ivi,vr∈H1(D),vi∈H1(D)

If we now pose in(18) v = Ā = A^r –iAⁱ, we immediately have

(19)k2||Ar||L2(D)2+||Ai||L2(D)2+∫∂DgrAr+giAids+i∫∂DgiAr−grAids=||∇Ar||L2(D)2+||∇Ai||L2(D)2+iβk∫∂D(Ar)2+(Ai)2ds,ork2||A||L2(D)2+∫∂DgA¯ds=||∇A||L2(D)2+iβk||A||L2(D)2.

We remind that the computational domain we assume in this section is rectangle, and our discretization mesh consists of squares with the edge of size h. We have m finite volumes along one direction of the domain D, and n along the another direction, so we have n×m finite volumes. Thus we have

m(p)=h2,m(σpq)=h,dpq=h.

We can express our numerical scheme in a similar way as it was done for the continuous equation. We obtain for p ∊ T_h,int

(20)k2uph2+∑q∈N(p)(uq−up)=0,

and for p ∊ T_h,ext

(21)k2uph2+∑q∈N(p)int(uq−up)+∑qext∈N(p)ext(uqext−up)=0.

We can approximate the boundary condition by

(22)uqext−uph+ikβupq=gpq,

where

upq=uqext+up2.

Substituting u_qext – u_p in (21) from (22) we have

k2uph2+∑q∈N(p)int(uq−up)+∑qext∈N(p)extgpq−ikβupqh=0.

Now we multiply each of the equation (20) and (21) by ū_p = upr−iupi, and add them together. Then for the second term using the usual finite volume property [6] we have

k2∑p∈Thup2h2+∑σpq∈Eextgpqu¯ph=∑σpq∈Eintuq−up2+ikβ∑σpq∈Eextupqu¯ph.

which can be rearranged in the form

k2∑p∈Thup2h2+∑σpq∈Eextgpqu¯pqh=∑σpq∈Eintuq−up2+ikβ∑σpq∈Eextupq2h+I,

where

I=∑σpq∈Eextgpq−ikβupq(u¯pq−u¯p)h.

First we notice

(u¯pq−u¯p)=12(u¯qext−u¯p),

and using boundary approximation (22) for real and imaginary part we have

(u¯pq−u¯p)=12g¯pqh+ikβhu¯pq.

Thus for the term I we obtain

I=12h∑σpq∈Eextgpq−ikβupqg¯pqh+ikβhu¯pq=12h∑σpq∈Eextgpq2+k2β2upq2+2kβgpqrupqr−gpqiupqih.

Now we use the discrete H¹ seminorm defined in [6]

(24)||u||1Th=∑σ∈Eintm(σ)dσ(Dσu)212,

where D_ςu = |u_p-u_q|,ς = σ_pq. If we now use the definition of a constant numerical solution (8), we can write

(25)k2||uh||L2(D)2+∑σpq∈Eextgpqu¯pqh=||uhr||1Th2+||uhi||1Th2+ikβ∑σpq∈Eextupq2h+h2∑σpq∈Eextgpq2+k2β2upq2+2kβgpqrupqr−gpqiupqih.

4 Numerical Scheme for Regular Mesh

We use the same computational domain as in the previous section. Now we derive numerical schemes for the real and imaginary part of the complex valued approximation function.

We have two types of the linear algebraic equations. The first belongs to the finite volumes from the set T_h,int, and second for the finite volumes from the set T_h,ext. In the first case we have for both, real and imaginary part of the numerical solution, the following equations

(26)k2uprh2+∑q∈N(p)(uqr−upr)=0,k2upih2+∑q∈N(p)(uqi−upi)=0.

For the second case we have some sides (one or two) of the finite volume that belongs to the boundary of the domain D. Here we use the boundary condition and the approximation described above

(27)uqextr−uprh−kβuqexti+upi2=gpqr,uqexti−upih+kβuqextr+upr2=gpqi.

From these equations we can easily eliminate the unknown values uqextr,uqexti:

uqextr=4hgpqr+2h2kβgpqi+4hkβupi+(4−h2k2β2)upr4+h2k2β2,uqexti=4hgpqi−2h2kβgpqr+(4−h2k2β2)upi−4hβkupr4+h2k2.

Substituting these values into (15) for p ∈ T_h,ext we have

(28)k2uprh2+∑q∈N(p)int(uqr−upr)+∑qext∈N(p)extAupr−Bupi+Cgpqr+Dgpqi=0,k2upih2+∑q∈N(p)int(uqi−upi)+∑qext∈N(p)extAupi+Bupr+Cgpqi−Dgpqr=0,

where

(29)A=−4+k2h2β24+k2β2h2,B=4kβ2h4+k2β2h2,C=4h4+k2β2h2,D=2kβ2h24+k2β2h2.

In this way we obtain linear system of algebraic equation with unknowns upr,upi,p∈Th.

Figure 4

Exact solution (left real part, right imaginary part) for the wavenumber 25rad/m

5 Numerical experiments

This section describes the results of the code, which solves the Helmholtz equation by the Finite volume method on the square domain of size 1 metre.

Figure 7

Numerical solution (left real part, right imaginary part) for the wavenumber 25rad/m, n = 60

5.1 Experiment 1 - boundary conditions with β = 1

First to be presented is the solution of (3)-(5) with β = 1 in (5), so Robin boundary conditions are prescribed for all sides of domain. For the beginning we must use the exact solution used in [3]

(30)u(x,y)=ei(k1x+k2y)=cos⁡(k1x+k2y)+isin⁡(k1x+k2y),

where the values k₁, k₂ are given by

(31)k1=kcos⁡θ,k2=ksin⁡θ.

The source function g was set so as to comply with the exact solution (30). Firstly we present results for θ=π2 , for two different wavenumbers. Figure 1 shows the exact solution for k = 10rad/m, extra for real and imaginary part. Figure 2 depicts the numerical solution to the Helmholtz equation for the number of discretizing points n = 10 (i.e. 100 finite volumes), and the Figure 3 for n = 40. It is clear that with the finer discretization the results are more accurate.

Figure 1

Exact solution (left real part, right imaginary part) for the wavenumber 10rad/m

Figure 2

Numerical solution (left real part, right imaginary part) for the wavenumber 10rad/m, n = 10

Figure 3

Numerical solution (left real part, right imaginary part) for the wavenumber 10rad/m, n = 40

Next figures are dedicated to the bigger wavenumber. It is known that higher frequencies (i.e. bigger wavenumbers) require finer discretization, if we want to get results close to the exact solution [7]. Figure 4 shows the exact solution for the wavenumber 25rad/m. Figures 5 and 6 depict the numerical solution for 10 and 40 discretizing points.

Figure 5

Numerical solution (left real part, right imaginary part) for the wavenumber 25rad/m, n= 10

Figure 6

Numerical solution (left real part, right imaginary part) for the wavenumber 25rad/m, n = 40

It is clear that it is more difficult to approximate this function. Figure 7 shows plots for n = 60, where the results get very close to the exact solution.

The L2 error was calculated by the formula

(32)∑((upr−uexactr)2m(p)2+(upi−uexacti)2m(p)2).

Here upr and upi are the numerically calculated values. uexactr and uexacti are the precise values calculated from the exact solution. Table 1 shows the values. The experimental order of convergence was calculated using the formula

(33)L2errorh<Chα,

Table 1

Values of the L2 error, θ=π2

	k= 10rad/m			k=25rad/m
n	10	40	60	10	40	60
L2	0.2653	0.0143	0.0063	1.3518	0.2324	0.1000
error

where h is the length of the finite volume. α will be then calculated by

(34)α=log2⁡L2errorhL2error2h

and from the theory it is known that its value is expected to be converging to 2. Table 2 shows the results for the wavenumber 10rad/m.

Table 2

Values of EOC for the wavenumber 10rad/m, θ=π2

n	L2 error	α
10	0.265255	2.17192
20	0.058864	2.03935
40	0.014320	2.01011
80	0.003555

Next figures 8 – 14 are for the value of θ=π4. First we show the exact solution and the results of the numerical experiments for the wavenumber 10rad/m ( Figure 8 - 10). Lastly, figures 11 - 14 depict the solutions for the wavenumber 25rad/m.

Figure 8

Numerical solution (left real part, right imaginary part) for the wavenumber 10rad/m

Figure 9

Numerical solution (left real part, right imaginary part) for the wavenumber 10rad/m, n = 10

Figure 10

Numerical solution (left real part, right imaginary part) for the wavenumber 10rad/m, n = 40

Figure 11

Exact solution (left real part, right imaginary part) for the wavenumber 25rad/m

Figure 12

Numerical solution (left real part, right imaginary part) for the wavenumber 25rad/m, n = 10

Figure 13

Numerical solution (left real part, right imaginary part) for the wavenumber 25rad/m, n = 40

Figure 14

Numerical solution (left real part, right imaginary part) for the wavenumber 25rad/m, n = 60

Presented figures and tables (3 and 4) show, that the behaviour is very similar to the previous value of θ.

Table 3

Values of the L2 error, θ=π4

	k= 10rad/m			k=25rad/m
n	10	40	60	10	40	60
L2	0.1425	0.0082	0.0036	1.5407	0.1238	0.0541
error

Table 4

Values of EOC for the wavenumber 10rad/m, θ=π4

n	L2 error	α
10	0.142520	2.09584
20	0.033340	2.02356
40	0.008200	2.00635
80	0.002041

5.2 Experiment 2 - mixed boundary conditions

The final part of the paper shows the results of program with changed boundary conditions. The domain is divided in two parts

(35)∂D=∂DRobin∪∂DNeumann.

The boundary conditions on three sides of the boundary ∂D_Robin are prescribed as in the previous case, where β = 1 in (5). For part of the right side of the domain

(36)∂DNeumann=∂D∖∂DRobin,

(37)∂DNeumann={[x,y]∈∂Ω;x=1;0.25<y,0.75},

the value of β = 0 in (5), so zero Neumann boundary conditions were obtained

(38)∂A∂n=0.

This combination of Robin and Neumann boundary conditions can represent for example a column (i.e. a hard wall barrier) standing in the free space. Exact solution is the same as in the previous case (30). The following tables 5 and 6 show the values of L2 error and EOC.

Table 5

>Values of the L2 error for case with changed boundary conditions, θ=π2

	k= 10rad/m			k=25rad/m
n	10	40	60	10	40	60
L2	0.2757	0.0150	0.0066	1.2876	0.2376	0.1024
error

Table 6

Values of EOC for the wavenumber 10rad/m for the case with changed boundary conditions, θ=π2

n	L2 error	α
10	0.275654	2.15734
20	0.061793	2.04152
40	0.015010	2.01061
80	0.003725

6 Conclusion

We have studied the numerical solution to the Helmholtz complex-valued equation. Our numerical solution is obtained using classical Finite volme method. For discretization of the boundary condition, which is of Robin type, we have used aditional exterior finite volumes by eliminating values on them. Properties of the weak solution and numerical solution are derived. Numerical experiments of various cases with changing boundary conditions show experimental order of convergence for the numerical solution to the exact one.

Acknowledgement

This work was supported by VEGA 1/0728/15.

References

[1] Chandler-Wild S., Langdon S., Boundary element methods for acoustics, Lecture notes, University of Reading, Department of Mathematics, 2007Search in Google Scholar

[2] Handlovičová A., Riečanová I., Roozen N. B., Rigid piston simulations of acoustic space based on Finite volume method, Proceedings of Aplimat, 15-th Conference on Applied Mathematics, Institute of Mathematics and Physics, Faculty of Mechanical Engineering STU in Bratislava, 2016, ISBN 978-80-227-4531-4Search in Google Scholar

[3] Babuška I. M., Ihlenburg F., Paik E. T., Sauter S. A., A Generalized Finite Element Method for Solving the Helmholtz Equation in two Dimensions with Minimal Pollution, Comput. Methods Appl. Mech. Engrg. 128,1995, 325-35910.1016/0045-7825(95)00890-XSearch in Google Scholar

[4] Riečanová I., Finite volume method scheme for the solution of Helmholtz equation, Proceedings of Aplimat, 14-th Conference on Applied Mathematics, Institute of Mathematics and Physics, Faculty of Mechanical Engineering STU in Bratislava, 2015, ISBN 978-80-227-4314-3Search in Google Scholar

[5] Sauter S. A., Schwab C., Boundary Element Methods, Springer-Verlag, Berlin, 2011, 56110.1007/978-3-540-68093-2Search in Google Scholar

[6] Eymard R., Gallouet T., Herbin R., Finite Volume Methods, Handbook of Numerical Analysis, Elsevier, 2007Search in Google Scholar

[7] Riečanová I., Solution of the Helmholtz equation by the Boundary element method and the Finite volume method, Proceedings of Advances in Architectural, Civil and Environmental Engineering, faculty of civil Engineering, Slovak University of Technology in Bratislava, 2014, 29-36, ISBN 978-80-227-4301-3Search in Google Scholar

Received: 2016-6-8

Accepted: 2016-10-11

Published Online: 2016-12-12

Published in Print: 2016-1-1

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.