Boundary conditions

In solving the fullwave model, the boundary conditions at an interface of different materials play an important role, and will significantly determine the solutions of a set of nonlinear partial equations. Hence we should assume carefully appropriate boundary conditions in order to obtain correct solutions.

In practical calculations, it is often essential to introduce artificial boundaries to limit the area of computation. Here we develop a systematic method for obtaining a hierarchy of local boundary conditions at these artificial boundaries. These boundary

conditions not only guarantee stable difference approximations, but also minimize the (unphysical) artificial reflections that occur at the boundaries.

At the edges of the grid, the stencil applied at points in the interior of the grid typically cannot be used. Some type of rule for handling the edge points is required.

There are several different types of boundary conditions:

(1) Dirichlet: , where xi is a boundary point. Sometimes instead of

vanishing of f, the unknown may be equal to some given function at the boundary.

|_xi 0 f =

|_xi f n

∂ =

(2) Neumann: ∂ 0, where n represents the coordinate that is normal to the boundary. Implementing this boundary condition on the grid using a forward difference approximation for the derivative leads to the relationship x₁= , x₂ where x is a point on the grid boundary. ₁

(3) Mixed: A linear combination of the function f and its normal derivative are set to a constant.

(4) Absorbing boundary condition (ABC): This type of boundary condition is very important in applications of the finite difference method, and is adopted by us in our simulation.

Most electromagnetic problems involve unbounded regions, which cannot be modeled computationally. One option is to use one of the above boundary

conditions and make the simulation region very large, and terminate the simulation before reflections from the boundary perturb the solution in the region of interest. The drawback of this approach is that the larger the simulation region, the greater the computational cost of the simulation. A better approach is to use a boundary condition that absorbs waves and reflects as little energy as possible.

This is the computational analogue of an anechoic chamber. There are several types of ABCs, including:

♦ One-way wave equation. These are easy to implement but imperfect in 2D and higher dimensions.

♦ Perfectly matched layer with loss.

♦ Surface integral equation on the boundary (MOM-TDIE).

Boundary conditions may also be required at material interfaces inside the simulation region.

Hyperbolic PDEs

Here we consider the simplest hyperbolic PDE - the 1D wave equation for function f:

( ) ( )

2 2

2 2 2

, 1 ,

f x t f x t

x c t

∂ ∂

∂ = ∂ (3.7.1)

One physical problem that is modeled by this PDE is a planar time-domain current

source that varies in intensity only in the x direction. Since the source varies only in the x direction, the radiated electric field also only varies in the x direction, and the components and both satisfy Eq. (3.7.1) where the source term is zero. In order to apply the finite difference method to the wave equation, difference approximations for the derivatives are required. The stencil for the second derivative in x is

where the t dependence of u(x; t) is suppressed for brevity. The stencil for the right-hand side of (3.2) is very similar. Substituting difference approximations into the

wave equation leads to

( ) ( ) ( ) ( ) ( ) (

The difference equation becomes

2 1

fm⁺ to obtain an explicit finite difference method:

Now, we can solve this for

+ ⎡ ⎤ −

This algorithm is known as the finite difference time domain method (FDTD), because a time coordinate is involved.

Because the wave equation involves time, part of the boundary condition required for the finite difference approach is actually a boundary condition in time, or an initial condition. Since the PDE involves a second order time derivative, initial conditions at two time steps are required. This means that f_m¹ and f_m² must be specified as given functions. One common situation is the initial condition f_m¹ = f_m²

= 0, and a source is applied at one of the spatial boundaries of the region.

( )

f t = ^{and f G t}

( )

A simplest boundary condition is Dirichlet: ,

where [0, G] is the simulation region. The Neumann condition is implemented by setting the endpoint value equal to the point next to it: f₂ⁿ = f₁ⁿ, and similarly for the right boundary point.

An absorbing boundary condition can be obtained by discretizing the one-way wave equation at the endpoints of the region. This is known as the Mur boundary condition.

At the right-hand side, we enforce the PDE

( )

( )

c

f x t f x t

x t

∂ ∂

∂ = − ∂ (3.7.6)

( )

(

)

f x t = f x− t

This equation has solutions of the form , which is a wave of arbitrary shape moving to the right as time increases. This allows waves to move out of the simulation region without reflection. We have to be careful in discretizing this

equation, because the approximations for the spatial and time derivatives in the one-way wave equation need to be evaluated at the same point.

There are two types of sources that can be used in the FDTD method, hard sources and soft sources. A hard source simply sets the value of the field at one or more grid points equal to a specified function of time, and so is a type of Dirichlet boundary condition. This corresponds to an EM problem in which the electric field at some point is known, and we wish to find the values of the radiated field at other points. One property of a hard source is that waves propagating towards the source are reflected by the source. A soft source corresponds to an impressed electric current. In order to allow for a soft source, we must rederive (3.7.2) from Maxwell’s equations. If we take the curl of Faraday’s law and substitute in Ampere’s law, we obtain

2

where the light speed , E and J are the electric field and current density,

respectively. If we use the vector calculus identity

( )

E E

−∇×∇× + ∇ ∇i = ∇ E

i i

(3.7.8)

and assume that the permittivity is constant and the net electric charge is zero, so that , then we arrive at the wave equation

0

Provided that there is a net charge density ρ in space, it will be obtained that

/

If we consider two dimensional case (i.e. x and z direction), then Eq. (3.7.10) can be

written as

For a problem in which the current density vector is in one direction only, and the source varies also only in one direction, this reduces to a 1D wave equation of the same form as Eq. (3.7.1) but with a forcing function determined by the current source.