The integral equation method for electromagnetic scattering problem at oblique incidence

Applied Numerical Mathematics 62 (2012) 860–873

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

The integral equation method for electromagnetic scattering problem at oblique incidence Haibing Wang a,b,∗ , Gen Nakamura a a b

Department of Mathematics, Hokkaido University, Sapporo, 060-0810, Japan Department of Mathematics, Southeast University, Nanjing, 210096, PR China

a r t i c l e

i n f o

Article history: Received 15 September 2011 Received in revised form 24 January 2012 Accepted 25 February 2012 Available online 13 March 2012 Keywords: Electromagnetic scattering Oblique incidence Solvability Integral equations method Numerics

a b s t r a c t We consider the scattering of electromagnetic waves scattered by an infinitely long impedance cylinder at oblique incidence, which is modeled as a system of a pair of the two-dimensional Helmholtz equations with coupled oblique boundary conditions. The solvability of such a scattering problem is proven by using the boundary integral equation method. By expressing the scattered fields in the form of single-layer potentials, our oblique scattering problem is transformed to a system of two integral equations. It is not a usual Fredholm system of the second kind as that in the case of normal incidence, since the system involves the tangential derivatives of the single-layer potential. By relating it to the Cauchy integral operator, we show that this system of operators is of Fredholm type with index 0. Therefore, the solvability of the integral system follows from the uniqueness of its solutions due to the Fredholm theory. A numerical scheme for solving the integral equations is also presented with some numerics. The numerical results illustrate the validity and efficiency of the proposed method. © 2012 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Scattering theory has been an active research field in engineering and applied mathematics. The electromagnetic scattering process is generally described by a boundary value problem for Maxwell’s equations defined in an infinite domain. Both the theoretical analysis and numerical methods for such a problem have been extensively studied, see [3,4,9,10]. Of particular interest in scattering theory is the scattering of a time harmonic electromagnetic wave by an infinitely long cylinder. In this case, the mathematical formulation of the scattering process is reduced to a two-dimensional problem. In particular, if the cylinder is illuminated by a TM or TE polarized plane wave at normal incidence, the corresponding electromagnetic scattering can be reformulated as a boundary value problem for the two-dimensional Helmholtz equation [2]. However, when the incident plane wave is obliquely incident at the cylinder parallel to z axis, the mathematical formulation of the scattering process consists of two boundary value problems for a pair of the two-dimensional Helmholtz equations in R2 . They are for the z components of the electric and magnetic fields, and the components are coupled through the boundary conditions even if the incident wave is a pure TM or TE polarized plane wave. Many researchers investigated the numerical methods for solving this kind of oblique scattering problems, see, for example, [1,5,11,14–18]. When the cross-section of the cylinder is of special shape, it is natural to seek the analytic solution by using the technique of separation of variables. For example, if the cross-section of the cylinder is circular, one can express the z components of scattered electric field and magnetic field as infinite series of the Hankel functions, and then determine the coefficients from the boundary conditions. For the cylinder with arbitrary cross-section shape, a moment method and a

*

Corresponding author at: Department of Mathematics, Hokkaido University, Sapporo, 060-0810, Japan. E-mail addresses: [email protected] (H. Wang), [email protected] (G. Nakamura).

0168-9274/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2012.02.006

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

861

finite element method are designed in [1] and [11], respectively. Nevertheless these numerical methods are verified to be efficient and accurate, the solvability of the forward scattering under considerations has not been addressed too much. This paper is devoted to the theoretical analysis and numerical study of the oblique scattering problem by using the integral equation method. More explicitly, we consider the scattering of a time harmonic electromagnetic wave by a cylinder at oblique incidence. We assume that the cylinder is not perfectly conducting but that does not allow the electromagnetic wave to penetrate deeply into the cylinder, which leads to a so-called Leontovich impedance boundary condition. Such a scattering problem is governed by a system of two Helmholtz equations with coupled boundary conditions. It should be emphasized that the boundary conditions involve tangential derivatives of the scattered fields. By expressing the scattered fields as the single-layer potentials, the scattering problem is transformed to two boundary integral equations. Note that the tangential derivative of single-layer is not a compact operator in usual function spaces, the solvability of resulting integral system cannot be proven in a simple way as that in the case of normal incidence. In the present work, we show that the resulting integral operator A of the boundary integral equations is Fredholm with index Ind( A ) = 0. We first decompose A as A = B + P , where P is compact and B involves Cauchy singular integrals. It turns out that B is Fredholm of index zero. Then, it follows from the Fredholm theory that A is also Fredholm and its index is zero. Hence, for establishing the solvability of the integral system, we only need to justify the uniqueness of its solutions. To generate stable numerical solutions of the oblique scattering problem, we propose a numerical scheme for solving the resulting integral formulation. Noting that the kernels of integral operators are singular, we decompose the kernels and extract their singularities in the form of Hilbert integral operator and Symm’s integral operator. Using trigonometric interpolation, the singular integrals can be efficiently approximated. Since our integral equation method is applicable to the obliquely incident scattering problem by the cylinder with cross-sections having arbitrary shape, it is expected that our present work can be widely used in some practical situations. The rest of this paper is organized as follows. In Section 2, we show the well-posedness of the forward scattering problem under considerations, by using the integral equation method. The resulting system of integral equations is proven to be Fredholm with index 0. A numerical scheme for solving this system is proposed in Section 3. We present some numerical results in Section 4, showing the validity and efficiency of our method. Finally we give some conclusions in Section 5. 2. Well-posedness of direct scattering problem In this section, we are ready to establish the well-posedness of the direct electromagnetic scattering problem from an impedance cylinder at oblique incidence. We assume that the background medium is homogeneous and dielectric. Let



   Ei , Hi = ei (x, y ), hi (x, y ) e −i ωt −i β z

be the time harmonic incident plane electromagnetic wave, where ω > 0 is the frequency and β is a constant depending on the incident direction. The total field is denoted as (E, H). By expressing the x components and y components of both the electric field E and magnetic field H by their z components (e , h), we derive from the Maxwell’s equations for (E, H) and the Leontovich impedance boundary condition (ν × E) × ν = λ(ν × H) that the scattering problem under consideration is modeled as a boundary value problem for a pair of the two-dimensional Helmholtz equations. More explicitly, let D be the cross-section of the cylinder, and (e s , h s ) be the z components of the scattered waves, where the time dependence and z dependence have been suppressed. Then we have the following system

e s + k2 e s = 0 in R2 \ D , s

2 s

(2.1)

2

h + k h = 0 in R \ D , s

2

k s ∂e λβ ∂ h +i e − = ∂ν ω ω ∂ τ β ∂ es ∂ hs λk2 s +i h + = ∂ν μω μω ∂ τ   √ ∂ es − ike s = 0, lim r r →∞ ∂r   √ ∂ hs − ikh s = 0, lim r r →∞ ∂r

λ

(2.2)

s

f1

on ∂ D ,

(2.3)

f2

on ∂ D ,

(2.4)

r = |x|,

(2.5)

r = |x|,

(2.6)

where ν and τ are the unit outward normal vector and tangent vector to ∂ D, respectively. Here  and μ are positive constants, which denote the permittivity and the permeability of the background medium, respectively. We assume that the impedance λ(x) for x ∈ ∂ D is positive. √ Let θ ∈ (0, π ) be the incident angle with respect to the negative z axis. Then we have β = ω μ cos θ and k =



√

μω2 − β 2 = ω μ sin θ . We remark that f 1 and f 2 can be given in terms of the incident wave (e i , hi ) as follows:

862

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

 k2 i ∂ ei λβ ∂ h i , +i e − ∂ν ω ω ∂ τ  i  ∂h β ∂ ei λk2 i f2 = − +i h + . ∂ν μω μω ∂ τ 

f1 = − λ

(2.7) (2.8)

We now begin with the uniqueness of solutions to the system (2.1)–(2.6). Theorem 2.1. If the impedance satisfies λ > 0, then the problem (2.1)–(2.6) admits at most one solution. Proof. It suffices to show that if e s and h s solve (2.1)–(2.6) for f 1 = f 2 = 0, then e s = h s = 0 in R2 \ D. Let Ωr be a disk with radius r centered at the origin and contain D in its interior. Set D r = Ωr \ D. If f 1 = f 2 = 0, then Green’s identity and the impedance boundary conditions (2.3), (2.4) imply

 e

s ∂e

s

∂ν

 ds =

e

s ∂e

∂ν

∂D

∂Ωr

e

∂D

ik

=

2

s

 ds +

 s 2

∇ e + e s e s dx

Dr





=

s

i



k

2

λω

es



 s 2 β ∂ hs ∇ e − k2 e s 2 dx + ds + ω ∂ τ

2 β (1/λ) e s ds +

ω

Dr

 e

ω

∂D

s ∂h

s

∂τ

 ds +



 s 2 ∇ e − k2 e s 2 dx

(2.9)

Dr

∂D

and

 hs

∂ hs ds = ∂ν

 hs

∂D

∂Ωr

∂ hs ds + ∂ν 



hs i

= ∂D

=i

k

2



μω

 Dr

λk

2

μω

hs −

 s 2

∇ h + h s h s dx 

 s 2 β ∂ es ∇ h − k2 h s 2 dx + μω ∂ τ

2 β λ h s ds −



hs

μω

∂D

Dr

∂ es ds + ∂τ





 s 2 ∇ h − k2 h s 2 dx.

(2.10)

Dr

∂D

If we note that

 es

 ∂ hs ∂ es ds = − h s ds, ∂τ ∂τ

∂D

∂D

then we see from (2.9) and (2.10) that

     2 2

∂ es ∂ hs k2  (1/λ) e s + λ h s ds  0,   es ds + μ hs ds = ∂ν ∂ν ω ∂Ωr

∂Ωr

(2.11)

∂D

where a denotes the imaginary part of a complex number a. On the other hand, from the radiation condition it follows that

2    s   s 2 s ∂e ∂e s + k2 e s 2 + 2k e s ∂ e ds = ∂ν ∂ ν − ike ds → 0, ∂ν ∂Ωr

∂Ωr

2    s   s 2 s ∂h ∂h s + k2 h s 2 + 2k h s ∂ h ds = ∂ν ∂ ν − ikh ds → 0, ∂ν

∂Ωr

∂Ωr

as r → ∞, which imply

lim

r →∞

s 2      s 2  s s ∂e ∂h 2 s ∂e s ∂h + μk2 h s 2 ds = −2k lim    +  k2 e s + μ e ds + μ h ds . r →∞ ∂ν ∂ν ∂ν ∂ν

∂Ωr

∂Ωr

∂Ωr

(2.12)

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

863

Hence, we have from (2.11) and (2.12) that

 lim

r →∞



s e ds = 0,

s h ds = 0.

lim

r →∞

∂Ωr

(2.13)

∂Ωr

Then, Rellich’s lemma immediately yields e s = 0 and h s = 0 in R2 \ D. Thus the proof is complete.

2

We next prove the existence of the solution to (2.1)–(2.6) by the integral equation method. To this end, we define the fundamental solution to the Helmholtz equation in R2 with wave number k by

Φ(x, y ) :=

(1 ) 

i 4



k|x − y | ,

H0

x, y ∈ R2 , x = y ,

(2.14)

(1)

where H 0 is the Hankel function of the first kind of order zero. Then, we introduce the following boundary integral operators defined by

 S [ϕ ](x) := 2

Φ(x, y )ϕ ( y ) ds( y ),

x ∈ ∂ D,

∂D



K [ϕ ](x) := 2

∂Φ(x, y ) ϕ ( y ) ds( y ), ∂ ν (x)

x ∈ ∂ D,

∂Φ(x, y ) ϕ ( y ) ds( y ), ∂ τ (x)

x ∈ ∂ D.

∂D





H [ϕ ](x) := 2 ∂D

Some properties of these operators are given as follows, which play a key role in the proof of our main result in this section. Lemma 2.2. Let ∂ D be of class C 2,α . Then the operators S , K : L 2 (∂ D ) → L 2 (∂ D ) are compact, and H : L 2 (∂ D ) → L 2 (∂ D ) is bounded. Proof. By Theorem 3.6 of [4], the operators S , K : L 2 (∂ D ) → H 1 (∂ D ) are bounded. Hence S and K are compact from L 2 (∂ D ) into itself, and H is bounded from L 2 (∂ D ) into itself due to the boundedness of the tangential operator ∂τ : H 1 (∂ D ) → L 2 (∂ D ). We further define the operators

 K 0 [ϕ ](x) := 2

∂Φ0 (x, y ) ϕ ( y ) ds( y ), ∂ ν ( y)

x ∈ ∂ D,

∂Φ0 (x, y ) ϕ ( y ) ds( y ), ∂ τ ( y)

x ∈ ∂ D,

∂Φ0 (x, y ) ϕ ( y ) ds( y ), ∂ ν (x)

x ∈ ∂ D,

∂Φ0 (x, y ) ϕ ( y ) ds( y ), ∂ τ (x)

x ∈ ∂ D,

∂D



H 0 [ϕ ](x) := 2 ∂D





K 0 [ϕ ](x) := 2 ∂D

H 0 [ϕ ](x) := 2



∂D

where

Φ0 (x, y ) =

1 2π

ln

1

|x − y |

,

x, y ∈ R2 , x = y

is the fundamental solution to the Laplace equation. We have the following important property [13]: Lemma 2.3. The operators H 0 and K 0 have the relation:



K 0

2

 2 − H 0 = I ,

where I is the identity operator.

(2.15)

864

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

Proof. Let V be the Cauchy integral operator defined by

V [ϕ ]( z) :=



1

ϕ (t ) t−z

πi ∂D

z ∈ ∂ D.

dt ,

Then, it follows from Section 7.5 of [7] that

V = −K0 + i H0.

(2.16)

Combining Theorem 7.10 with Theorem 7.26 of [7], we have

V2 = I

in L 2 (∂ D ).

As a consequence, it holds that

K 02 − H 02 = I , and hence



K 0

2

 2 − H 0 = I ,

because K 0 and H 0 are adjoint operators of K 0 and H 0 , respectively.

2

We are now in a position to state the main result of this section: Theorem 2.4. If −k2 is not a Dirichlet eigenvalue of the Laplacian in D and λ ∈ C (∂ D ) is positive, then the system (2.1)–(2.6) has a unique solution (e s , h s ). Proof. We seek the solution to (2.1)–(2.6) in the form of single-layer potential



e s (x) =

Φ(x, y )ϕ1 ( y ) ds( y ),

x ∈ R2 \ D ,

(2.17)

Φ(x, y )ϕ2 ( y ) ds( y ),

x ∈ R2 \ D

(2.18)

∂D



s

h (x) = ∂D

ϕ1 , ϕ2 ∈ L 2 (∂ D ). Using the standard jump relations of single-layer potential, we derive from (2.3)

with unknown densities and (2.4) that

−ϕ1 + K [ϕ1 ] + i β

μω

k2

λω

S [ϕ1 ] −

H [ϕ1 ] − ϕ2 + K [ϕ2 ] + i

β

ω

λk2

μω

H [ϕ2 ] =

A :=

2

− I + K + i λkω S β

μω H



β

f 1,

S [ϕ2 ] = 2 f 2 .

Let A be the operator defined by

2

λ

ω H

(2.19) (2.20)



2

λk I − K − i μω S

.

Then, in terms of

 T

ϕ := (ϕ1 , −ϕ2 ) ,

f :=

2

λ

T f 1, 2 f 2

,

the system (2.19)–(2.20) becomes

A [ϕ ] = f ,

(2.21)

where the superscript T denotes the transposition. Now we proceed to prove the solvability of (2.21) in two steps. In the first step, we show that A is a Fredholm operator with index 0. To this end, we decompose A as

A = sin θ B + P with

(2.22)

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873



1

B :=

sin θ

μω



β

−I β

ω H 0

H

865

 (2.23)

I

0

and

2



ω ( H − H 0 )

K + i λkω S

P :=

β

2



μω ( H − H 0 )

λk − K − i μω S

β

 (2.24)

.

Since the kernels of H and H 0 have the same singularity, the operator

H − H 0 : L 2 (∂ D ) → L 2 (∂ D ) is compact, and hence

P : L 2 (∂ D ) × L 2 (∂ D ) → L 2 (∂ D ) × L 2 (∂ D ) is compact. We next show that B is of Fredholm type with index zero. Using Lemma 2.3, we derive that

B =

1

2

β2



I + μω2 ( H 0 )2

0

0

I + μω2 ( H 0 )2

sin2 θ

β2

 =

I 0

0 I

 +

1 sin2 θ

β2

μω2

( K 0 )2

0

0 β2

μω2

( K 0 )2

 .

Due to the compactness of K 0 , the operator B is Fredholm, and by the index theorem [8] we have that





Ind B 2 = 2 Ind( B ) = 0. It follows from (2.22) that A is also a Fredholm operator with

Ind( A ) = Ind( B ) = 0. In the second step, we justify the uniqueness of the solutions to (2.21). Then the solvability of (2.21) follows from the Fredholm alternative theorem. Let



ϕ0 := ϕ10 , −ϕ20

T

be the solution to the homogeneous form of integral equation (2.21), that is, f = 0. Then the potentials e 0s , h0s with densities ϕ10 , ϕ20 solve the homogeneous form of (2.1)–(2.6). By Lemma 2.1, we have

e 0s = 0,

h0s = 0 in R2 \ D ,

(2.25)

h0s = 0 on ∂ D .

(2.26)

and hence

e 0s = 0,

In addition, the potentials e 0s , h0s defined by (2.17) and (2.18) also satisfy the Helmholtz equation in D. Since −k2 is not a Dirichlet eigenvalue of the Laplacian in D, we have from (2.26) that

e 0s = 0,

h0s = 0 in D .

(2.27)

Using the jump relation of single-layer potential, we have



∂ e 0s ∂ν

+

 −

∂ e 0s ∂ν

−

 = −ϕ10 ,

∂ h0s ∂ν

+

 −

∂ h0s ∂ν

−

= −ϕ20 ,

(2.28)

where the indices + and − denote the standard limits from R2 \ D and D onto the boundary ∂ D, respectively. Combining (2.25), (2.27) and (2.28) yields

ϕ10 = ϕ20 = 0 on ∂ D . The proof is complete.

2

(2.29)

866

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

3. Numerical scheme for solving the integral system In this section, we develop an effective numerical scheme for solving the integral equations (2.21). Since (2.21) involves singular integrals, we cannot directly apply the usual quadrature rules to discretize it. To generate a stable numerical solution, we decompose the kernels using periodic extractions of the singularities. Note that the operators S and K have weakly singular kernels, while H has a Cauchy singular kernel. We are led to deal with the periodic Hilbert integral operator and Symm’s integral operator [12]. A quadrature method based on trigonometric interpolation is used to generate a stable approximation. We now describe the decomposition scheme of the kernels. Assume that the closed curve ∂ D is smooth enough, and has the parametrization

    ∂ D := x := x(r ) = x1 (r ), x2 (r ) : r ∈ [0, 2π ]

(3.1)

with 2π -periodic functions xi (r ), i = 1, 2. We define

 2 2 x (r ) := x (r ) + x (r )

  ˜ r ) := ψ(x) = ψ x(r ) , ψ(

1

2

for a given density ψ and x = x(r ) ∈ ∂ D. By direct calculations, we have the parametric expressions of K , S and H :





K [ψ] x(t ) =



2π

1

|x (t )|



(3.2)

0

2π





˜ r ) x (r ) dr , K(t , r )ψ(



˜ r ) x (r ) dr , S(t , r )ψ(

S [ψ] x(t ) =

(3.3)

0

  H [ψ] x(t ) =

2π

1

|x (t )|





˜ r ) x (r ) dr H(t , r )ψ(

(3.4)

0

with the kernel functions

K(t , r ) := S(t , r ) := H(t , r ) :=

ik 2 i 2

(1 ) 

H0

ik 2

H 1(1) (k|x(r ) − x(t )|)



n(t ) · x(r ) − x(t )

|x(r ) − x(t )|

,







H 1(1) (k|x(r ) − x(t )|)

k x(r ) − x(t ) ,

p (t ) · x(r ) − x(t )

|x(r ) − x(t )|

,

where











n(t ) = x 2 (t ), −x 1 (t ) = x (t ) ν x(t ) ,











p (t ) = x 1 (t ), x 2 (t ) = x (t ) τ x(t ) .

Here we have made use of ( H 0 ) = − H 1 with H 1 being the Hankel function of the first kind of order one. It is evident that S(t , r ) has the logarithmic singularity at t = r, and H(t , r ) involves the Cauchy-type singularity. K(t , r ) is finite, but its derivatives are singular at t = r. In terms of the singularities of kernels, we decompose the functions K(t , r ), S(t , r ) and H(t , r ) as follows: (1)

(1)

 K(t , r ) = K1 (t , r ) ln 4 sin

2

 S(t , r ) = S1 (t , r ) ln 4 sin2

2 t −r

 H(t , r ) = H1 (t , r ) ln 4 sin2 with

t −r

 + K2 (t , r ),



2 t −r 2

(1)

+ S2 (t , r ),  +

1 2π

cot

r −t 2

+ H2 (t , r )

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873



r )−x(t )|) − 2kπ n(t ) · [x(r ) − x(t )] J 1 (|kx(|xr ()− , t = r , x(t )| t = r,

(3.5)

− 21π J 0 (k|x(r ) − x(t )|), t = r , t = r, − 21π ,

(3.6)

r )−x(t )|) , t = r , − 2kπ p (t ) · [x(r ) − x(t )] J 1 (|kx(|xr ()− x(t )| 0, t = r,

(3.7)

K1 (t , r ) =

 S1 (t , r ) =



867

0,

H1 (t , r ) =

where J 0 and J 1 are the Bessel functions of order zero and order one, respectively. The kernels K2 , S2 and H2 are defined by

 K2 (t , r ) = K(t , r ) − K1 (t , r ) ln 4 sin2

 S2 (t , r ) = S(t , r ) − S1 (t , r ) ln 4 sin

2

t −r 2

t −r

 ,

2

 H2 (t , r ) = H(t , r ) − H1 (t , r ) ln 4 sin2

 ,

t −r

 −

2

1 2π

cot

r −t 2

with diagonal terms

1 x

(t ) · n(t )

K2 (t , t ) =

|x (t )|2

2π

S2 (t , t ) =

,

i

−

2

γ 1 k|x (t )| − ln , π π 2

H2 (t , t ) = −

1 x

(t ) · p (t ) 2π

|x (t )|2

,

where γ is the Euler’s constant. According to the above decompositions, we are faced with the following two singular integrals:

2π cot

r −t 2

˜ r ) dr and ψ(

0

2π  ln 4 sin2

t −r 2



˜ r ) dr . ψ(

0

Using trigonometric polynomials for approximation [6,14], we have

2π cot

r −t 2

˜ r ) dr ≈ ψ(

2N −1 

(N )

Tj

˜ t j ), (t )ψ(

tj =

j =0

0

2π  ln 4 sin2

t −r 2



˜ r ) dr ≈ ψ(

2N −1 

˜ t j ), R j (t )ψ( (N )

jπ N

,

tj =

j =0

0

jπ N

with the weights (N )

Tj

(t ) =

π N

tj −t 1 − (−1) cos( Nt ) cot , j

2

(N )

R j (t ) = −

2π N

 N −1  1 m =1

m

cos m(t − t j ) +

1 2N

 cos N (t − t j ) .

Other regular integrals involved in (3.2)–(3.4) possess periodic smooth kernels, and hence they can be approximated accurately by simple rectangular rule. For more details, we refer to [6,7]. 4. Numerical examples To illustrate the performance of the method described in this paper, we now present some numerical results and analyze the errors of the numerical solutions for two model problems. In the first example, we focus on a special case where the scattered waves e s and h s are expressed exactly by Φ(·, z). Such a model problem enable us to check the accuracy of the proposed method and demonstrate the efficient convergence property of the numerical solution. Noting that the incident plane wave is of great importance in practical situations, we consider the scattering problem from an obliquely incident plane wave in the second example. 4.1. A special case: The exact solution is known We proceed by constructing a model problem such that its true solution can be expressed analytically. For this purpose, we let z1 and z2 be two points inside D, and then define two boundary functions by

868

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

(1 ) 

f 1 = −λkH 1

 ν (x) · (x − z1 )

k|x − z1 |

|x − z1 |

+i

( 0) 

k2

ω



k|x − z1 | +

H1

λβ

ω

(1 ) 

 τ (x) · (x − z2 )

k|x − z2 |

kH 1

|x − z2 |

and

 ν (x) · (x − z2 )

(1 ) 

f 2 = −kH 1

k|x − z2 |

|x − z2 |

+i

λk2

μω



( 0) 

k|x − z2 | −

H1

(1)

β

μω

(1 ) 

kH 1

 τ (x) · (x − z1 )

k|x − z1 |

|x − z1 |

.

(1)

It is easy to verify that e s (x) = H 0 (k|x − z1 |) and h s (x) = H 0 (k|x − z2 |) satisfy the following system:

e s + k2 e s = 0

in R2 \ D ,

h s + k2 h s = 0

in R2 \ D ,

∂ es k2 s λβ ∂ h s +i e − = f 1 on ∂ D , ∂ν ω ω ∂ τ β ∂ es λk2 s ∂ hs +i h + = f 2 on ∂ D , ∂ν μω μω ∂ τ  s  √ ∂e lim r − ike s = 0 r = |x|, r →∞ ∂r   √ ∂ hs lim r − ikh s = 0 r = |x|. r →∞ ∂r λ

According to the asymptotic behavior of the Hankel function, the far-filed patterns of e s and h s can be expressed as

−4ie i π /4 −ikxˆ ·z1 e e ∞ (ˆx) = √ , 8π k

(4.1)

−4ie i π /4 −ikxˆ ·z2 e h∞ (ˆx) = √ 8π k

(4.2)

with the observation direction xˆ = (cos t , sin t ) ∈ S, t ∈ [0, 2π ], where S is the unit circle in R2 . On the other hand, we obtain from (2.17) and (2.18) that the far-field patterns of e s and h s can be expressed in the form



e i π /4

∞

e (ˆx) = √

8π k

e

−ik xˆ · y

e i π /4

ϕ1 ( y ) ds( y ) = √

∂D

8π k

2π









e −ikxˆ ·x(r ) ϕ˜ 1 (r ) x (r ) dr ,

(4.3)

0

and

e i π /4

∞

h (ˆx) = √ 8π k

 e

−ik xˆ · y

e i π /4

ϕ2 ( y ) ds( y ) = √

∂D

8π k

2π

e −ikxˆ ·x(r ) ϕ˜ 2 (r ) x (r ) dr ,

(4.4)

0

respectively, where ϕ˜ j (r ) := ϕ j ( y ) = ϕ j (x(r )) for y = x(r ) ∈ ∂ D. Hence, solving the integral system (2.21) by the numerical scheme presented in Section 3, we can compute the far-field patterns of e ∞ and h∞ via (4.3) and (4.4), respectively. And then, we compare them with the exact solutions (4.1) and (4.2). Let D be a kite-shaped domain with the boundary expression

x1 (t ) = 2.0 cos t + 1.3 cos 2t − 1.3,

x2 (t ) = 2.5 sin t ,

t ∈ [0, 2π ].

The geometric shape is described in Fig. 1. z1 = (0.2, 0.2) and z2 = (−0.3, −0.3) are two points inside D which are used to construct our model problem. In the following, we set the constants  , μ and ω equal to 1. The boundary impedance is taken as λ = 2. In numerical implementations, we choose the equidistant collocation points by setting

tj =

π N

j,

j = 0, . . . , 2N − 1.

We give the numerical solutions of e ∞ and h∞ for N = 32 in Figs. 2 and 3, respectively. The point-wise relative errors of |e ∞ | and |h∞ | are also displayed in Fig. 4. It is clear that our numerical solutions approximate the exact solutions very well, with a relative error order 10−8 . This fact greatly confirms the effectiveness of the proposed method. To demonstrate the convergence behavior, we provide the absolute errors of |e ∞ | and |h∞ | for different N at four special points (see Table 1). As N increases, the accuracy of the numerical solutions is remarkably improved.

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

869

Fig. 1. The geometry of D. z1 and z2 are two points inside D

Fig. 2. The comparison of numerical solution with exact solution for e ∞ (solid line: exact solution computed via (4.1); open circles: numerical solution generated by (4.3), N = 32).

Table 1 Absolute error of |e ∞ | for different N. t

N =8

N = 16

N = 32

0

3.967E-003 7.419E-003 1.502E-003 5.051E-003

3.364E-005 6.309E-006 1.729E-005 2.634E-006

1.204E-008 2.770E-008 9.880E-009 3.840E-008

π /2 π 3π /2

Table 2 Absolute error of |h∞ | for different N. t

N =8

N = 16

N = 32

0

6.648E-003 3.114E-003 6.484E-003 5.425E-003

1.181E-005 9.601E-006 1.903E-005 4.952E-006

6.723E-009 1.780E-008 3.992E-009 1.961E-008

π /2 π 3π /2

870

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

Fig. 3. The comparison of numerical solution with exact solution for h∞ (solid line: exact solution computed via (4.2); open circles: numerical solution generated by (4.4), N = 32).

Fig. 4. The point-wise relative errors of |e ∞ | and |h∞ | for N = 32.

4.2. A general case: Incident plane wave Since the incident plane waves are frequently encountered in real situations and the corresponding far-field patterns are the input data of inverse scattering problems, we now turn our attention to the general electromagnetic scattering problems from an impedance cylinder by obliquely incident plane waves. Let

E i (x; d, p ) = √

1

 k20

1 curl curl pe ik0 x·d = √ (d × p ) × de ik0 x·d ,

(4.5)



1 1 H i (x; d, p ) = √ curl pe ik0 x·d = √ d × pe ik0 x·d i μk0 μ

(4.6)

√

be the incident electric field and magnetic field with the wave number k0 := ω μ , where the time dependence has been suppressed. In particular, we let the incident wave be a TM polarized plane wave at an angle θ with respect to the negative

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

871

Fig. 5. The norm of scattered wave e s in a rectangular domain [−5, 5] × [−5, 5]. Table 3 Numerical solutions of |e ∞ | for different N. t

N = 16

N = 32

N = 64

0

(−0.0578917509, −0.2311300062) (−1.087921356, 1.000149236) (0.0369829469, −0.1329154063) (−0.3548502594, 0.2120746645)

(−0.0578692277, −0.2311266321) (−1.087919982, 1.000145114) (0.0369834195, −0.1328905717) (−0.3548758326, 0.2121061668)

(−0.0578692882, −0.2311266321) (−1.087919977, 1.000145112) (0.0369833970 − 0.1328906008) (−0.3548757903, 0.2121060956)

π /2 π 3π /2

z axis, and the plane of incidence has an angle ξ with the x axis. Then, the incident direction and the polarization can be expressed as

d = (sin θ cos ξ, sin θ sin ξ, − cos θ) and

p = (cos θ cos ξ, cos θ sin ξ, sin θ), respectively. From (4.5) and (4.6), it follows that



cos θ cos ξ E (x; d, p ) = √ cos θ sin ξ e ik0 sin θ (x cos ξ + y sin ξ ) e −izk0 cos θ ,  sin θ

 sin ξ 1 i H (x; d, p ) = √ − cos ξ e ik0 sin θ (x cos ξ + y sin ξ ) e −izk0 cos θ . μ 0 1

i

Hence, for the TM polarized plane wave, the fields e i and h i in (2.7) and (2.8) are given by

sin θ e i (x, y ) = √ e ik0 sin θ (x cos ξ + y sin ξ ) ,

(4.7)

h i (x, y ) = 0.

(4.8)



Here we note that e i (x) corresponds to a plane wave in R2 with wave number k = k0 sin θ and incident direction d˜ = (cos ξ, sin ξ ) ∈ S. Let the domain D and the constants  , μ, ω be the same as those in last subsection. We next report some numerical results for the above TM polarized incident plane wave. In Figs. 5 and 6, we display the distributions of the norm of the scattered fields |e s | and |h s | for ξ = π /2, respectively. Since the exact solution is unknown, we cannot show the error of the numerical solution. But the convergence property can be easily observed from Tables 3 and 4, where the numerical solutions for different numbers N are given. 5. Conclusions In this paper, we investigated the scattering of a time harmonic electromagnetic wave by an impedance cylinder at oblique incidence. The mathematical formulation of such a scattering problem consists of two Helmholtz equations with

872

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

Fig. 6. The norm of scattered wave h s in a rectangular domain [−5, 5] × [−5, 5]. Table 4 Numerical solutions of |h∞ | for different N. t

N = 16

N = 32

N = 64

0

(0.0176771085, 0.2108770319) (−0.0052048371, −0.0155859986) (−0.0485263019, 0.0149741799) (−0.0202981068, −0.0261858698)

(0.0176690387, 0.2108695535) (−0.0051730401, −0.0156084996) (−0.0484980448, 0.0149762552) (−0.0203014856, −0.0261725008)

(0.0176690575, 0.2108695665) (−0.0051731200, −0.0156084590) (−0.0484980989, 0.0149762446) (−0.0203014754, −0.0261725285)

π /2 π 3π /2

coupled boundary conditions for the z components of electric and magnetic fields. The principal difficulties for both theoretical analysis and its numerical realization lie in the fact that tangential derivatives are involved in the boundary conditions. Using the integral equation method, we proved the solvability of the oblique scattering problem under consideration. The key is to verify that the system of integral equations associated with the scattering problem via single layer potentials is Fredholm with index 0. To generate a stable numerical solution, we decompose the kernels and extract the singularities. By this, we are led to deal with the Hilbert integral operator and Symm’s integral operator, which can be efficiently approximated by using trigonometric polynomials. In addition, we presented two numerical examples to show the effectiveness of the proposed scheme. By constructing a model problem such that its exact solution can be expressed analytically, the accuracy and convergence behavior of the numerical solution are greatly illustrated. Acknowledgements The authors would like to thank the referees for their careful reading and valuable suggestions. Especially, we highly appreciate one referee who gave us a way to simplify the proof of Theorem 2.4 and suggested us using the numerical scheme presented in Section 3. This scheme turns out to be more effective than that used in our original manuscript, and makes the paper much improved. This work is supported by Grant-in-Aid for Scientific Research (B) (22340023) of the Japan Society for the Promotion of Science. The author H.B. Wang is also partially supported by NSF of China (Nos. 11071039, 11161130002). References [1] A.C. Cangellaris, R. Lee, Finite element analysis of electromagnetic scattering from inhomogeneous cylinders at oblique incidence, IEEE Transactions on Antennas and Propagation 39 (1991) 645–650. [2] F. Cakoni, D. Colton, Qualitative Methods in Inverse Scattering Theory, Springer-Verlag, Berlin, 2006. [3] D. Colton, R. Kress, Integral Equation Methods in Scattering Theory, John Wiley & Sons, New York, 1983. [4] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer-Verlag, Berlin, 1998. [5] M. Lucido, G. Panariello, F. Schettiho, Scattering by polygonal cross-section dielectric cylinders at oblique incidence, IEEE Transactions on Antennas and Propagation 58 (2010) 540–551. [6] R. Kress, On the numerical solution of a hypersingular integral equation in scattering theory, J. Comput. Appl. Math. 61 (1995) 345–360. [7] R. Kress, Linear Integral Equations, Springer-Verlag, New York, 1999. [8] W. McLean, Strong Elliptic Systems and Boundary Integral Equations, Cambridge University Press, London, 2000. [9] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, Oxford, 2003. [10] J.C. Nédélec, Acoustic and Electromagnetic Equations, Springer-Verlag, New York, 2001. [11] R.G. Rojas, Scattering by an inhomogeneous dielectric/ferrite cylinder of arbitrary cross-section shape-oblique incidence case, IEEE Transactions on Antennas and Propagation 36 (1988) 238–246.

H. Wang, G. Nakamura / Applied Numerical Mathematics 62 (2012) 860–873

873

[12] J. Saranen, G. Vainikko, Periodic Integral and Pseudodifferential Equations with Numerical Approximation, Springer-Verlag, Berlin, 2002. [13] G. Schmidt, Integral equations for conical diffraction by coated gratings, J. Int. Eqns. Appl. 23 (2011) 71–112. [14] J.L. Tsalamengas, Exponentially converging Nyström methods applied to the integral-integrodifferential equations of oblique scattering/hybrid wave propagation in presence of composite dielectric cylinders of arbitrary cross section, IEEE Transactions on Antennas and Propagation 55 (2007) 3239– 3250. [15] N.L. Tsitas, E.G. Alivizatos, H.T. Anastassiu, D.I. Kaklamani, Optimization of the method of auxiliary sources (MAS) for oblique incidence scattering by an infinite dielectric cylinder, Electr. Eng. 89 (2007) 353–361. [16] J.R. Wait, Scattering of a plane wave from a circular dielectric cylinder at oblique incidence, Canadian J. Phys. 33 (1955) 189–195. [17] J. Yan, R.K. Gordon, A.A. Kishk, Electromagnetic scattering from impedance elliptic cylinders using finite difference method, Electromagnetics 15 (1995) 157–173. [18] H.A. Yousif, A.Z. Elsherbeni, Oblique incidence scattering from two eccentric cylinders, J. Electromagnetic Waves Appl. 11 (1997) 1273–1288.