An FDTD model for low and high lightning generated electromagnetic fields

PDF

Download

Open document

Flip pages

Contact Us

Copy asset link

Request this asset

Transcript (if available)

Content

An FDTD Model for Low and High Lightning Generated Electromagnetic Fields

By

Svetlana Knyazeva

_______________________________________________________________________

A Thesis Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
MASTER OF ARTS
(APPLIED MATHEMATICS)

May 2013

Copyright 2013 Svetlana Knyazeva

1

Table of Contents

Abstract ..................................................................................................................................................................... 2
I. Introduction ...................................................................................................................................................... 3
II. Literature survey ............................................................................................................................................... 4
III. Physics of the Problem ...................................................................................................................................... 6
Physical Geometry ................................................................................................................................................ 6
Field Equations ...................................................................................................................................................... 7
Maxwell’s equations in Cylindrical Do main .......................................................................................................... 8
Boundary Condition .............................................................................................................................................. 9
Ionosphere assumptions ..................................................................................................................................... 10
Cold Plasma Parameters ..................................................................................................................................... 10
IV. Numerical Domain .......................................................................................................................................... 11
2-Dimensional Computation Grid ....................................................................................................................... 11
Discrete Equations .............................................................................................................................................. 11
Finite-Difference Time-Domain Method ............................................................................................................. 12
Iteration Coefficients .......................................................................................................................................... 14
Boundary Conditions at PML and on the Earth’s Surface ................................................................................... 15
Air-Ground Interface ........................................................................................................................................... 16
Absorbing Boundary Condition – Nearly Perfect Match Layer ........................................................................... 16
Simulation Parameters and Inputs...................................................................................................................... 17
V. Numerical Validation ...................................................................................................................................... 19
VI. Conclusion ....................................................................................................................................................... 21
VII. References ....................................................................................................................................................... 22
Applendix A ............................................................................................................................................................. 24
Iteration Matrices ............................................................................................................................................... 24
Appendix B............................................................................................................................................................... 29
Matlab Code Organization .................................................................................................................................. 29
Appendix C ............................................................................................................................................................... 31
Electrons and Ions Densities ............................................................................................................................... 31

2

Abstract
One of the natural sources of the electric current is the lightning discharge in the ionosphere. The
electro-magnetic wave produced by intense lightning, like sprites, can significantly alter electrostatic,
electromagnetic, acoustic, radar, and radio-frequency measurements in ionosphere and on the ground.
There are numerous studies performed on the sprites and other types of lightning discharge activities.
This study concentrates on very low frequency (VLF, 3-20 kHz) and extremely low frequency (ELF, 3-3000
Hz) electromagnetic waves radiated by lightning currents. These frequencies have very low attenuation
rate and make the ground a good enough conductor where the radio wave arrives at the surface
entering the ground nearly independent of the incident angle. The Earth curvature can be neglected for
frequencies up to ~12 kHz since the wave incident angles are very small. Method being used is Finite-
Difference Time-Domain (FDTD) in the two-dimensional cylindrical simulation environment. Surface
impedance boundary condition is used for the air-ground interface; and Nearly Perfect Match Layer is
implemented as absorbing boundary condition on the top and sides of the simulation domain.

3

I. Introduction
The Earth’s ionosphere which is the layers of the atmosphere at altitudes between 50 to 2,000
kilometers above the Earth surface plays an important role in wireless communication. This is due to the
fact that, as its name suggests, the ionosphere contains substantial concentration of ionized gas. The
physical processes of ionization of neutral gas by the solar radiation and the reverse process of
recombination following sun set are well understood for the dominant region in terms of ion density,
the F-region, at altitudes of 120 to 2,000 kilometers of the ionosphere. Substantial efforts have been
successful in developing Numerical Weather Prediction (NWP) capabilities for the monitoring and the
forecasting of the F-region ionosphere [15], [14]. However, the layers of the ionosphere that often have
the most impact on the wireless communication on Earth are the E and D-regions of the ionosphere.
These regions of the ionosphere are located at altitudes of 50 to 100 kilometers above the Earth’s
surface. The densities of ions are substantially lower in the E and D-region than in the F-region. The low
ion density presents significant challenges in the observation of the ion densities in these regions.
In the measurement of the ion contents of any layer of the ionosphere, the primary technique consists
of transmitting radio signals between two locations and using the distortion such as phase delay of the
signal to infer the electron content along the path of transmission of the signal. The most common
measurements of the ionosphere are the signal delays from received Global Positioning Satellites (GPS)
by either ground based or space based receivers. However, the E and D-region of the ionosphere have
little effect on the GPS signal at over 1GHz frequency. As the Appleton-Hartree model of refractivity of
magnetized plasma indicates [16] that the frequencies of radio wave most effective at measuring the
electron content of the E and D regions of the ionosphere are from 50-1,000 Hz for Extremely Low
Frequency (ELF) and 3 to 30 kHz for Very Low Frequency (VLF). There is very limited number of
transmitters and receivers operating at the VLF and ELF range dedicated for the monitoring of the
ionosphere. It is therefore useful to make use of naturally occurring lightning events as the source of the
radio wave for the measurement of the E and D-region ionosphere. In fact, electromagnetic waves
radiated by these lightning currents propagate between the ground and ionosphere with low
attenuation rates. The long-range detectability of these signals enables ionosphere remote sensing, and
remote lightning discharge current waveform measurement [1].
The low frequency of the VLF and ELF radio waves also translates into very long wavelengths. The
commonly used ray-tracing technique for modeling radio wave propagation through plasma may be
inaccurate. This is due in part to the fact that the ray-tracing equation relies on the Eikonal
approximation of the solution of the Maxwell’s equation and requires the “small wavelength”
assumption. Alternatively, it is possible to model the space between the ionosphere and the surface of
the Earth as a wave guide. We can then solve the Maxwell equation to determine the direction and the
strength of the magnetic and electric fields perturbation due to the propagation of the lightning signal.
The primary objective of this thesis is to examine the mathematical techniques for modeling the
disturbance of the electric and magnetic fields due to lightning. The work presented in this thesis follows
closely the modeling approach developed by S.A. Cummer and his co-authors. In particular, we use the
idealized cylindrical domain and model for the lightning discharge of current in Cummer’s work. The
focus of our effort is on the following areas:
4

1. Understand all basic assumptions needed for the modeling of the disturbance of electric and
magnetic fields created by a lightning strike;
2. Study the derivation and the properties of the numerical schemes used to approximate the
solution of the Maxwell’s equation;
3. Fine-tune the parameters used in the Perfectly Matching Layers (PML) approach for the
modeling of an infinite domain for the wave propagation, and
4. Develop a prototype simulation packages to evaluate the accuracy of the approximate solutions
and the sensitivity of the solutions to different model parameters.
In Chapter II, we present a brief review of the literature in the modeling electro-magnetic waves created
by the lightning. In Chapter III, an overview of the Maxwell’s equation and the relevant models for
current in a magnetized plasma medium is provided. In Chapter III and IV we shall present the numerical
scheme used in our model, as well as, the boundary conditions including the PML approach for our
model. We shall present the numerical validation results in Chapter V which is followed by the
conclusions.
II. Literature survey
The big portion of literature relating to this work was written by Dr Steven A. Cummer at the Duke
University and his students (people.ee.duke.edu/~cummer/). While he was a student himself, he
studied different approaches to electromagnetic fields modeling in the ionosphere. He summarized
these studies in his 1997 paper “An Analysis of New and Existing FDTD Methods for Isotropic Cold
Plasma and a Method for Improving Their Accuracy” [3]. Cummer demonstrates a first interest in the
very low and extremely low frequencies electromagnetic waves in “Ionospheric D region remote
sensing using VLF radio atmospherics” paper in 1998 [4].
This paper is based on the publication by Wenyi Hu and Steven A. Cummer “An FDTD Model for Low
and High Altitude Lightning-Generated EM Fields” [1]. They validate their model by comparing results to
the mode theory results.
The main simulation algorithm – Finite-Difference Time domain was first introduced by Kane Yee in 1966
and published in IEEE Transactions on Antennas and Propagation. Yee describes the basis of the FDTD
numerical technique for solving Maxwell’s curl equations directly in the time domain on a 3D space grid.
Since first introduction FDTD was improved, adopted for other coordinate systems and modified to work
with a range of stability criteria and assumptions. One of the Cummer’s papers [3 ] explores the use of
FDTD and demonstrates good results that are compatible with popular mode theory for electromagnetic
waves.
In the last few years, the finite-difference time-domain (FDTD) method has been extended to modeling
electromagnetic wave propagation and interactions with cold plasmas [3], [10]. Many FDTD models
treat ionosphere as a simple nonmagnetized isotropic medium. However, for accurate investigations,
effects introduced by the Earth’s magnetic field on the ionospheric plasma cannot b e ignored. Thus, a
magnetized cold plasma ionospheric medium must be used in the FDTD simulations. The ionosphere
can be regarded as inhomogeneous cold plasma with Earth’s magnetic field superposed as long as the
5

energy generated by lightning is not high enough to modify the medium. Although the electric fields
produced by intense lightning discharges can modify the ionosphere, this sort of effect is localized, and
the assumption of a linear cold plasma medium should not affect the simulation results significantly [1],
[3].
A few two-dimensional (2D) and three-dimensional (3D) FDTD models have been published that include
the magnetic effect on the ionosphere. S.A. Cummer proposed a 2D cylindrical-coordinate FDTD model
to study electromagnetic wave propagation in the Earth-ionosphere waveguide [10], and then extended
it to explore lightning-generated electromagnetic wave behavior in ionosphere and to test the sprite
initiation theory [1]. The stability condition for his model is E-J collocation method, unlike the H-J
collocation method used in other models. This method collocates E (electric field) and J (current density)
in time and space. E-J collocation method is independent of medium properties and remains the same
for free space [1], [3]. This is a very important characteristic of Cummer’s FDTD algorithm, and this
paper uses the same stability condition.
Source current is the lightning discharge and the formula is taken from the following source: A.S. Dennis
and E.T. Pierce, “The return stroke of the lightning flash to Earth as a source of VLF atmospherics,” Radio
Sci., vol. 68D, p. 777, 1964.
The simulation results are very depended on the parameters of the ionosphere like ion densities, plasma
frequencies, gyrofrequencies, and collision frequencies. These parameters are provided by the Ion
Reference Ionosphere (http://iri.gsfc.nasa.gov/), and empirical standard model of the ionosphere, based
on all available data sources.
Since no one has computing resources to simulate full domain of the wave propagation there is a need
of boundary conditions. The goal of a boundary condition is to absorb all outward traveling waves
without producing any reflections. Although this goal is not entirely realistic since all boundary
conditions produce some degree of back reflection, the better the boundary condition the lower the
reflections will be. The nearly perfectly matched layer (NPML), a versatile absolute boundary condition
(ABC) that is simple to implement in complicated medium was adopted in this work.
The development of highly absorbing boundary conditions (ABC) has been one of the most active areas
in FDTD research over the last few years. Kane Yee started this work when he published FDTD algorithm
[7]. As stated earlier, the goal of an ABC is to absorb all outward traveling waves without producing any
reflections. The ABCs can be implemented using the Bayliss-Turkel annihilation operator, the Engquist-
Majda one-way wave equations, Mur boundary conditions, Liao's extrapolation, and the perfectly
matched layer (PML) method [2]. The PML method was originally formulated by Berenger in 1994 for
use with Maxwell's equations, represents one of the most significant advances in FDTD development,
since its conception by Kane Yee in 1966 [11]. Cummer explored perfectly matched layer and nearly
perfectly matched layer, compared them, and published results in “A Simple, Nearly Perf ectly Matched
Layer for General Electromagnetic Media” paper in 2003 [ 2] and “The Nearly Perfectly Matched Layer is
a Perfectly Matched Layer” paper in 2004 [9].

6

III. Physics of the Problem
Physical Geometry
Lightning in nature occurs in the atmosphere which is bounded from below by the Earth’s surface and
unbounded otherwise horizontally and from above. However the numerical simulation of the infinite
horizon problem in 3-dimensional space is computationally intensive. For the purpose of a proof of
concept study of radio wave propagation, we follow the approach used by Cummer to reduce the
problem to a 2-dimensional domain by considering a radially symmetric domain centered at the location
of the lightning discharge. In order to simulate the unbounded nature in both radial and vertical
directions, we use the approach of adding Perfectly Matched Layers (PML) to absorb the expending
disturbance wave from the lightning.
A schematic depiction of the modeling domain for our work is shown in Figure 1 below. A vertical
electric source (lightning) between ionosphere and the ground is located at the center. From the center
electromagnetic waves propagate outwards therefore cylindrical coordinate system is the most
convenient [1] to represent the components of electric and magnetic fields, as well as, the modeling
domain. We wish to calculate the electric and magnetic fields at the all locations between the source of
the lightning and the sampling point where disturbances of electrical and magnetic fields are measured
at the end of the modeling domain.

(a) (b)
Figure 1 (a) Cylindrical Computational Domain. PML stands for Perfectly Match Layer. (b) Slice of the
domain starting from the center (lightning discharge) to the end.
One main difference between 2D cylindrical model and fully 3D model is how Earth’s magnetic field is
treated. The Earth’s magnetic field has a significant impact on the radio wave behavior in the
ionosphere because the electric currents caused by the charged particles are affected by the
background magnetic field. In reality, the background magnetic field is not azimuthally symmetric and
thus cannot be simulated exactly in 2D cylindrical coordinates used in our FDTD code. However, a
reasonable approximation can be made in which the 1D background magnetic field variation along the

Lightning Discharge
Ground
PML
PM
PML
PML
Ground
Lightning Discharge
Sampling Point
Z
r
7

propagation path is rotated to provide azimuthal symmetry [1]. As shown in Figure 2 below, the radial
component

of the background geomagnetic field in the plane defined by the lightning source and the
sample point is replicated in all azimuthal directions. The same approximation is made for the and z
components of Earth’s magnetic field [1].

Figure 2. Diagram’s of Earth’s magnetic field radial component.
Field Equations
The electromagnetic fields in cold plasma are described by Maxwell’s equations coupled to equations
for current derived from the Lorentz equations of motion of the charged particles in the medium in
response to the wave electric field and an ambient static magnetic field from the Earth.

The fields in the medium are described by Maxwell’s equation s:

Differential Form Integral Form
Ampere’s Law

Faraday’s Law

where E = electric field, H = magnetic field, J = current density,

is vacuum permittivity defined by the
formula

and

is vacuum permeability. The integral form of these
equations is satisfied for any oriented surface and its oriented boundary .

Electric Currents follow the Lorentz’s equation of motion of the free electrons and ions in the medium in
response to the electric field and an ambient static Earth’s magnetic field :

Lorentz’s Equation
  E b J
q
q
J
t
J
pn E n Bn
n
n
n n
n 2
0
       



,

Source
Propagation Direction
Source
Propagation Direction
Side View Top View
8

where n = type of the particle (electron, positive ion, negative ion) and b
E
= unit vector in the direction
of the Earth’s magnetic field ,

= gyrofrequencies for each particle,

=collision frequencies,

=plasma frequencies.

A very common form of lightning is cloud-to-ground discharge (CG). This form is also the one that has a
well-defined current equation that is free of time stamp and spatial location; therefore it is easy to use
as a source current in this work [5]. We also assume almost perfectly vertical cloud-to-ground lightning
with uniform current discharge. Electric Current source is assumed to be uniformed along the lightning
discharge channel, located in the center of the modeling grid. Giving by:

Current Source     L t bt at
v
I t I J
source
/ ) exp( 1 ) exp( ) exp( ) (
0
0


       ,

where I
0
=20 kA, ν
0
= 8*10
7
m/s, γ = 3*10
4
1/s, a = 2*10
4
1/s, b = 2*10
5
1/s, and L = lightning discharge
channel length (meters) [5]. Amplitude of the source current is displayed in Figure 3 below. In our work,
we assume that all components of the source current density are identical, that is,

Figure 3. Lightning-strike electric current source at the center of the domain (discharge channel length is
1 kilometer).
Electric currents add up the following way:

Maxwell’s equations in Cylindrical Domain
The integral form of the Maxwell equations can be very useful in motivating for the formulation of the
finite difference approximation scheme. Consider a rectangular element in a cylindrical coordinate
system and 3 of the surfaces of element. Figure 4 illustrates the relationship between the total flux of
the magnetic fields through these surfaces and the line integral of the electric field along the boundaries
of these surfaces. Note that the red arrows represent the positive normal direction of the surfaces and
the corresponding direction of the boundary.
9

Figure 4. Application of Maxwell’s on full 3D cylindrical coordinate system . Vectors represent the
direction of the integral.
For an example, an approximate form of the Faraday’s law for the -component of the magnetic field
has the form

where

represents the average of the -component of the magnetic field over the surface area and

,

,

and

correspond to the average r-component and the z-component of the electric
field over the upper (+z), lower (-z), outer (+r) and inner (-r) sides of the boundary, respectively. Similar
approach can be made to approximate the Ampere’s Law. The integral form of these equations also
suggests that a nested grid should be used for the discretization of the electrical and magnetic fields as
we shall explain later.
Boundary Condition
The interpretation of the integral form of the Maxwell equation shown in Figure 4 can also be used to
determine the necessary boundary conditions for the well-posedness of the solutions. In fact since the
time derivatives of the electric and magnetic field through a surface element depend on the line integral
of the other field, the boundary values on the modeling domains are required for determining a solution.
On the outer, top and bottom of the modeling domain, Dirichlet boundary conditions are satisfied for all
components of the electric and magnetic fields. That is

and

at the top most and outer
most grid points of the PML region, as well as, on the Earth’s surface. Even though the notional center of
the cylindrical area is not an actual boundary, but we only solve the equation over one 2-dimension slice
of the cylinder, the center of the cylinder is a boundary of our modeling domain. We can find
appropriate boundary conditions using the geometric symmetry of the problem. Since we assume that
the lightning source is vertically uniform along the center of the cylinder and the fields are radially
symmetric with respect to the center, the following conditions are imposed at the center
Lightning Discharge
Lightning Discharge
Lightning Discharge
Orientation of the Surfaces and Arcs
10

The last set of conditions is also satisfied throughout the modeling domain.
Ionosphere assumptions
Based on the literature survey – the most interesting altitude range for calculating electro-magnetic
fields due to lightning is from ~55km to ~200km. In this case algorithm has to run through the whole
domain which is set from ground level to ~200km.
Due to the frequencies ranges given the effect of positive and negative ions is included. In many cases,
the effect caused by ions can be neglected because the mass of an ion is about 2000 to 6000 times
larger than an electron’s. Therefore at ELF/VLF, the field driven ion motion is usually much smalle r than
the electron motion. However, in the lowest ionosphere (below ~60km), the negative ion density is high
compared to electrons due to the complicated chemistry processes and dominates the conductivity of
the atmosphere at low altitudes [1]. The ions are O
+
, H
+
, He
+
, NO
+
, O
2
+
.
The effect of the wave magnetic field on the charged particle is neglected [1].
Cold Plasma Parameters
Due to the significant difference in the ionosphere parameters between nighttime and daytime this
simulation assumes nighttime profile. Nighttime profile and some physical aspects of the lightning
generation allow ionosphere to be treated as true cold plasma. Therefore, thermal motion is neglected
for the purpose of this simulation [1].
Plasma Frequency:

, where n=type of the particle (electron, positive ion, negative
ion),

is the charge of each type of particles,

is the mass of each species and

is the number
density of the particle.
Gyrofrequencies:

, where B is Earth’s magnetic field.
Collision frequencies:

(a) (b)
Figure 5. Ionosphere information taken from International Reference Ionosphere. (a) Electron and Ion
Density (the same values). (b) Collision frequency of electrons and ions.
11

IV. Numerical Domain
2-Dimensional Computation Grid
Our primary computational model is expressed in 2D azimuthally symmetric cylindrical coordinates.
Figure 6a shows the computational grid that is symmetric in directions. The bottom part is the ground,
which can be treated as a perfect electrical conductor (PEC) boundary for VLF frequencies. The
ionosphere is defined by number density of different charged particles and their collision frequencies;
electron and ion densities vary throughout the computational domain. The vertical and radial outer
regions are set as perfect matched layers as the absorbing boundaries [1]. The entire computational
domain is first mapped into discrete grids with desired spatial resolution. A vertical electric source is
located in the center of the cylindrical domain with adjustable length [5]. Then time iterations are to
solve the fields’ quantities at all the grids and all time steps.

(a) (b)

(c) (d)
Figure 6. (a) 2D cylindrical computational domain (b) grid layout
Discrete Equations
The task of the modeling becomes solving the governing equations everywhere in the computational
domain by approximating the curl operator as the finite discrete difference. For the calculations:

Ground PEC
PML Absorption oundary
PML
r
z
Vertical
Dipole
Source
Electron, ion densities and collision
frequencies vary throughout
the computational domain
Solve for E, H,
at each g rid point

r+1 r
z
z+1
r+1 r
z
z+1
r+1 r
z
z+1
12

0 



because all fields are assumed to be azimuthally symmetric. We first write down the differential
form of the Maxwell’s equation in component form.
Faraday’s Law
The differential form of the Faraday’s equation in cylindrical coo rdinate is given by

By separating the components, we obtain:
(1)

,
(2)

,
(3)

Ampere’s Law
Similarly, the differential form of the Ampere’s equation in cylindrical coordinate is given by

.
Component-wise forms of the above equation are given by

(4)

(5)

(6)

Equations (4)-(6) should be rewritten as E field update:
(4a)

(5a)

(6a)

Lorentz’s Equation
The Lorentz’s equation must be written for each of the 8 components of part icles in the ionosphere:
(7) – (15)   E b J
q
q
J
t
J
pn E n Bn
n
n
n n
n 2
0
       



,
where n represents the 3 species of particles in the ionosphere and x represents the 3 spatial directions.
Finite-Difference Time-Domain Method
The main simulation algorithm is the Finite-Difference time-domain (FDTD) method. The Electric field,
Magnetic field, and Electric current density are approximate at two sets of nested grid points respective.
13

As shown in Figures 6 (b)-(d), the electric field and current densities are defined at the 4 corner points of
the large yellow square and the magnetic field is defined at the 4 corner points of the large blue square.
A time marching scheme is used to propagate the states of the system in time [7]. For the calculations,
we recall that all partial derivative with respect to are zero because all fields are assumed to be
azimuthally symmetric.
The flavor of FDTD being used in this work is Direct Integration (DI) FDTD [3]. This method collocates E
(electric field) and J (current density) in time and space.

A time discrete version of the Faraday’s law is given by
(1b)

,
(2b)

,
(3b)

Similarly, a time discrete version of the Ampere’s Law FDTD equations using E -J collocation in time and
space is given by
(4b)

(5b)

(6b)

Using E-J collocation in time and space for Lorentz’s equation, we obtain

where

is a unit vector in the direction or Earth’s magnetic field.

FDTD version of the above Lorentz’s equation is broken down for 3 different type particles in 3
directions for 9 equations in total.

A time discrete version of the Lorentz’s Equations for the electron, n=e, in component form is given by

14

In the above equations, the terms involving solutions at time step N+1 are all placed to the left-hand
side of the equations. Therefore the scheme corresponds to a mixed explicit and implicit time marching
scheme. Equations for other particles are treated similarly. Detail expressions for these components are
given in the Appendix A.

Iteration Coefficients
Simulation deals with all three types of particles and 3 field components have to be calculated in 3
directions, which leave us with 15 variables in 15 equations. The whole system is divided into two
groups.

A(12x12), B(12x12) and C(12x6) are the coefficient matrixes that depend on the medium properties and
discretization parameters. For all equations above and below consider n=type of particle (electron,
positive ion, negative ion) and N= time step.

For the ease of computation in Matlab matrix A can be inverted and the equation is rewritten:

15

The entries in

and

are all the coefficients needed to implement the explicit FDTD
iteration and these coefficients are only needed to be calculated once for a time invariant system.
Please see details for A, B, and C matrices in the Appendix A. These coefficients vary with height so they
are precalculated based on the location and stored before main FDTD algorithm runs.
Equations (1), (2), and (3) are not used as iteration equations because E – J collocation in time and space
method is used. These equations are used to update H field and the Y vector at each grid cell.

Where

is the source current at time N+1/2 – it is calculated by given formula.
Boundary Conditions at PML and on the Earth’s Surface
Solutions to boundary conditions, both air-ground interface and absorbing boundary conditions, involve
Maxwell’s equations in the frequency domain. Through the inverse Fourier transform, general solutions
of Maxwell’s equations can be stated as linear combinations of single -frequency solutions.
Ampere’s Law

16

Faraday’s Law

Where is the input frequency of electromagnetic wave.
Air-Ground Interface
Applying FDTD method directly to the air-ground border is very computationally expensive since it’s
requires very small step size, especially for VLF waves. Approximation method used is the Surface
Impedance boundary condition, which takes the following assumptions into account: ground is a lossy
dielectric, which is important at the higher VLF frequencies and no fields inside the ground need to be
calculated. Method finds relationships between tangential electric fields and tangential magnetic fields
at the surface [12]. The frequency domain equations are:

Equations used in this simulation are simplified versions that assume the transmission angle of the wave
is always very close to

.








  ) ( ) (
'
t H
t
Z t E
r s 

0
0 '
/ ) (


 
i
i Z Z
s s
 
where

is the input frequency.
Absorbing Boundary Condition – Nearly Perfect Match Layer
In general the electromagnetic analysis of scattering structures often requires the solution of “open
region" problems, where the spatial domain of the computed fields is unbounded in one or more
coordinate directions. Obviously, one cannot compute the solution over such a space, due to limitations
in computer memory and computational time. As a result, it therefore becomes necessary to truncate
the computational domain in such a way as to make it appear infinite in extent. This is achieved by
enclosing the structure by a suitable output boundary that absorbs all outward traveling waves. Such a
boundary is referred to as an absorbing boundary, or an absorbing boundary condition (ABC), and its
function is to absorb all incident radiation and to suppress spurious back reflections of the outgoing
wave.
When modeling electromagnetic wave propagation there is a need of the boundary condition that
absorbs the wave without reflections and without abrupt cut-off at the border. In this model we are
dealing with cold plasma and absorbing boundary condition that is strong enough for such dispersive
and anisotropic medium is Nearly Perfect Match Layer (NPML). Another advantage of NPML is its
equations are very easy to integrate to different numerical methods including the numerical method in
this paper.

17

Where

and

are NPML conductivities in r and z direction and

.
Simulation Parameters and Inputs
Ions
In order to calculate plasma frequencies and gyrofrequencies we need mass and charge for each particle.
Mass and charge for the electron are well known, and therefore are easy inputs; it is not the same for
ions. There are 5 different ions that are present in the ionosphere (O
+
, H
+
, He
+
, NO
+
, O
2
+
) with two
different flavors: positive and negative [1]. We need average mass for ions and average positive and
negative charges. During the literature research for this work, no comprehensive information on masses
or charges of ions, or standard averages was found. For the simulations purposes the following
assumptions were used:
a) Since [1] states that mass of ion is 2000 to 6000 times larger than the mass of electron:
average mass of ion = 4000*mass of electron
b) negative ion charge=electron charge
c) positive ion charge=abs(electron charge)

Assumptions stated above can be changed in the code.
Electron and Ion density in ionosphere
For the plasma frequency calculations particle densities in the ionosphere are required. As stated earlier
this information is taken from International Reference Ionosphere organization (IRI). IRI provides
reference numbers based on the inputs such as location (latitude, longtitude), exact date (from 1958 to
2015), and time of the day.

IRI only explicitly provides electron densities in the ionosphere. Based on the literature survey and
images posted in Cummer’s papers, for the night time for the altitude above 80km electron and ions
densities are the same, and for below 80km, electron density is zero and ions density is the same as at
80km [1], [10]. Based on the inputs to IRI, the density at 80km was 7.182*10
7
m
-3
, the same density for
ions is assumed for the altitude below 80km.

For this research electron densities were extracted from the IRI based on the inputs specified in
Appendix C. The densities are stored in the Matlab function ParticalDensity.m. Figure 5(a) is the exact
graph of the densities.
Electron and Ion collision frequencies
Collision frequencies for each particle

are calculated by the following formulas [4]:

18

Where frequencies are the same for negative and positive ion and z is the altitude measured in
kilometers.
Earth’s magnetic field
For this purpose of this work Earth’s magnetic field is assumed to be the same in all three directions and
is equal to 5*10
-5
Tesla [1].

These magnitudes can be modified in the Matlab code as well.
Current Source
The lightning Current Source formula

is used in simulations [1],[5].
    L t bt at
v
I t I J
source
/ ) exp( 1 ) exp( ) exp( ) (
0
0


      
L is the discharge current length in meters. It is assumed to be 1km for all simulations; this values is
entered and can be changed in currentSource.m Matlab function.

FDTD algorithm was also run with a variation of the current source formula to study changes of the
electric currents away from the source.

Where c=3*10
8
(m/s) is the speed of light, ω is the input frequency and r is the radial distance from the
source.
Times
Cummer suggests to use maximum time step of 1.67*10
-6
seconds for a grid step size of 0.5 km or 1km
[1]. The total time to run the simulation depends on the input frequency. In order to calculate valid
electric and magnetic fields away from the current source algorithm should allow one full wave period.
Table 1 below states minimum times to run the simulation for some input frequencies.

Extremely Low Frequency Very Low Frequency
Frequency (Hz) Min Time (s)
3 0.333333333
100 0.01
300 0.003333333
1000 0.001
1500 0.000666667
3000 0.000333333

Frequency (Hz) Min Time (s)
3500 0.000285714
5000 0.0002
10000 0.0001
15000 6.66667E-05
20000 0.00005
30000 3.33333E-05

Table 1. Wave period=min time to run the simulation for the given frequencies.
Spatial dimension and step size
Spatial step size in both r and z directions are assumed to be 1km. While recommended size of Δz is
1/10 of the wave length and it should ideally vary with maximum altitude of the simulation[1], in this
work Δz is still assumed to be 1km for the ease of comparison the final results with these of Dr. St even
Cummer [1] and for the fact that our electron and ion densities are given for each kilometer of the
altitude.
19

Spatial step size can be changed in the Matlab code: FindIterationMatrices.m and mainFDTD.m.
Boundary Condition
There are two important parameters that are used in NPML boundary: conductivities in r and z
directions. The actual values for these parameters used in calculations were not stated by Cummer and
Hu, only theoretical background was described. Cummer in his 2003 paper provides suggestion for the
conductivities [2]:

Where r is the distance from the source. No suggestion for the top-layer of NPML was found. This
relationship was used as a starting point in the simulation for both side and top NPML.

V. Numerical Validation
The goal of this thesis is to reproduce detail equations for FDTD algorithm and all its inputs and to
calculate magnitudes of electromagnetic fields produced by lightning. This model is based upon the 2-D
cylindrical domain with magnetized cold plasma developed by Hu and Cummer [1]. Validation of the
accuracy of these calculations in Matlab code is very important. The simulations results of this paper are
compared to the simulation results produced by Hu and Cummer [1]. The exact results are not expected
because of slightly different inputs for ionosphere parameters and because exact dimensions of the
domain for the results were not always specified by Cummer [1]. However, results and their visual
representation are expected to be close.
Equations
Cummer and Hu stated general Maxwell’s equations and a compact form of their iteration algorithm
without deriving any final equation or the coefficient matrices [1]. This work represents the effort
completed to derive all 15 equations needed to implement FDTD simulation, together with integration
of equations for boundary conditions and cold plasma parameters. The final form of the iteration
algorithm derived here in Chapters 3 and 4 perfectly fits with the Cummer and Hu’s form [1].
Expected Values for Electric and Magnetics fields
To validate equations FDTD algorithm was run over the 100 time step of 1.67*10
-6
seconds each step for
frequency of 10kHz. Images below in Figure 7 are generated values for magnetic and electric fields at
different distances from the source at the altitude of 100 kilometers.

20

(a)

(b)

(c)

(d)

(e) (f)
Figure 7. (a)–(c) Magnetic Field in different directions (d)-(f) Electric Field in different directions.
21

Validation with FDTD simulation by Cummer.
Cummer and Hu have few plots in their paper that can be used as a visual reference of the calculated
values [1]. The luck of knowledge of exact inputs for the FDTD algorithm prohibitss exact comparison
between results.

VI. Conclusion
The prime outputs of this thesis are detailed equations used for FDTD algorithm and a Matlab code that
reproduces these equations and that can be reused for future work. The Matlab code functions and
inputs are described in the Appendix B. Code can be re-used, with inputs and reference numbers
modified.
There are several suggestions for the future work that came up during this research. One suggestion is
to move the location of the sample point to another place. It is currently located right at the boundary
layer (PML), where FDTD algorithm is being modified to generate infinity. Moving sample point towards
the end of the computational grid but not placing right at the border should theoretically produce better
results.
Iteration matrices A, B, and C all depend on inputs that are not clearly stated and use lots of
assumptions, like particles masses, densities, charges, and altitude step size. Calculating iteration
matrices based on the latitude-longitude ranges and comparing results from the same geographical
region would make it easier to validate any future electromagnetic calculations. Also, iteration matrices
as they are stated in this work tend to be ill-conditioned and should be either scaled up or scaled down,
depending on the inputs, to avoid stability and accuracy problems.
This work neglects Earth’s curvature, which seems to be numerically valid for extremely small frequency
ranges and for short distance propagation [1]. In general, Earth-flattening approximations would be a
good addition to this model and can be added to the Matlab code.
Despite several suggestions for the future additions, this work and the code can be used for the
electromagnetic fields calculations generated by lightning. It can also be used as a comparison or
validation tool for different approach on the same problem.
22

VII. References

[1] W. Hu, S.A. Cummer, “An FDTD Model for Low and High Altitude Lightning -Generated EM
Fields”, IEEE Transactions on Antennas and Propagation, Vol. 54, No5, May 2006.
[2] S.A. Cummer , “A Simple, Nearly Perfectly Matched Layer for General Electromagnetic Media”,
IEEE Microwave and Wireless Components Letters, Vol.13, No. 3, March 2003,
[3] S.A. Cummer , “ An Analysis of N ew and Existing FDTD Methods for Isotropic Cold Plasma and a
Method for Improving Their Accuracy”, IEEE Transactions on Antennas and Propagation, vol. 45, no. 3,
March 1997.
[4] S.A. Cummer, U.S. Inan, and T.F. ell, “Ionospheric D region remote sensing us ing VLF radio
atmospherics”, Radio Sci., vol. 33, num. 6, p. 1781, 1998.
[5] A.S. Dennis, E.T. Pierce, “The return stroke of the lightning flash to Earth as a source of VLF
atmospherics,” Radio Sci., vol. 68D, p. 777, 1964.
[6] J.P. erenger, “A perfectly matched layer for the absorption of electromagnetic waves,”
J.Comput.Phys., vol. 114, pp.185-200, 1994.
[7] K.S. Yee, “Numerical Solution of Initial oundary Value Problems Involving Maxwell’s Equations
in Isotropic Media, “ IEEE Transactions on Antennas and Propagation, vol.AP-14, p.302, 1966.
[8] F.L. Teixeira, W.C. Chew, “PML -FDTD in Cylindrical and Spherical Grids,” IEEE Microwave and
Guided Wave Letters, vol.7, num 9, 1997.
[9] W. Hu, S.A. Cummer, “The Nearly Peferctly Matched Layer is a Perfectly Matc hed Layer, “ IEEE
Antennas and Wireless Propagation Letters, vol. 3, 2004.
[10] S.A. Cummer, “Modeling Electromagnetic Propagation in the Earth -Ionosphere Waveguide,”
IEEE Transactions on Antennas and Propagation, vol.48, p.1420, 2000.
[11] . Schneider, “ Understanding the FDTD Method. Chapter 11 – Perfectly Matched Layer”,
www.eecs.wsu.edu/~schneidj/ufdtd/chap11.pdf
[12] S. armada, L.D. Rienzo, N. Ida, S. Yuferev, “The Use of Surface Impeda nce Boundary Conditions
in Time Domain Problems: Numerical and Experimantal Validation”, ACES Journal, vol. 19, no.2, 2004.
[13] L. Mandrake, . Wilson, C. Wang, G. Hajj, A. Mannucci, X. Pi, “A performance evaluation of the
operational Jet Propulsion Laboratory/University of Southen California Global Assimilation Ionospheric
Model PL/USC GAIM ”, Journal of Geophysical Research, vol. 110, A12306, 2005.
[14] G. Hajj, . Wilson, C. Wang, X. Pi, G. Rosen, “Data assimilation of ground GPS total electron
content into a physic-based ionospheric model by use of the Kalman filter”, Radio Science, vol. 39, 2004.
23

[15] C. Wang, G. Hajj, X. Pi, G. Rosen, . Wilson, “Development of the Global Assimilative Ionospheric
Model”, Radio Science, vol. 39, 2004.
[16] K. Davies, “Inospheric Radio”, Short Run Press Ltd., Exeter, UK 1990

24

Applendix A
Iteration Matrices
A MATRIX
E =Earth e = electron ip = ion positive in = ion negative

Plasma frequency (x = e, in, ip)

collision frequency (x = e, in, ip)

vector of the Earth’s magnetic field

25

A MATRIX Simplified
E =Earth e = electron ip = ion positive in = ion negative

Plasma frequency (x = e, in, ip)

collision frequency (x = e, in, ip)

vector of the Earth’s magnetic field

26

B MATRIX
E =Earth e = electron ip = ion positive in = ion negative

Plasma frequency (x = e, in, ip)

collision frequency (x = e, in, ip)

vector of the Earth’s magnetic fiel d

27

B MATRIX simplified
E =Earth e = electron ip = ion positive in = ion negative

Plasma frequency (x = e, in, ip)

collision frequency (x = e, in, ip)

vector of the Earth’s magnetic field

28

C MATRIX

29

Appendix B
Matlab Code Organization
Users should run the following two functions in the order stated:
1) FindIterationMatrices.m
Calculates A, B, C matrices, inverses of A and B matrices, A
-1
B and A
-1
C matrices.
Calls CollisionFrequency.m
PlasmaFrequency.m
Gyrofrequency.m
Inputs Outputs
Maximum time step A
-1
B matrix for each altitude step, saved in the text file and as a
global variable in Matlab domain
Magnitudes of Earth’s magnetic
fields in z, r, directions
A
-1
C matrix for each altitude step, saved in the text file and as a
global variable in Matlab domain
Δz = altitude step size
z = maximum altitude

2) mainFDTDcode.m
Runs the actual FDTD algorithm and saves electric and magnetic field values for each time step over the
whole computational domain.
For each time_step
For each r_step
For each z_step
Update X vector
Update H fields (Hr, Hz, Hphi)
Update Y vector – to be used in the next time step
Store all the X, Y, H values for each time_step.

Calls currentSource.m
Inputs Outputs
Maximum time step, total time of simulation Hr_matrix: Hr values in each grid point for each time
step
Maximum altitude z_max, Δz = altitude step size Hphi_matrix: Hphi values in each grid point for each
time step
Distance of the sample points =r_max, Δr = radial
step size
Hz_matrix: Hz values in each grid point for each time
step
Frequency X_matrix: X vector in each grid point for each time
step
A
-1
B and A
-1
C matrices that are defined globally. Y_matrix: Y vector in each grid point for each time
step
NPML conductivity power

30

The following functions are not run directly by the user but called by the either
FindIterationMatrices.m or by mainFDTDcode.m:
currentSource.m
CollisionFrequency.m
Collision frequencies

are calculated for each particle type.
PlasmaFrequency.m
Calls ParticleDensity.m
ParticleDensity.m
The table with particle densities extracted from IRI is stored here.
Gyrofrequencies.m

Miscellaneous functions used for graphing:
Graph_ParticalDensity.m
darpa_collisionf.m

31

Appendix C
Electrons and Ions Densities

International Reference Ionosphere organization (IRI) provides reference numbers based on the inputs
such as location (latitude, longtitude), exact date (from 1958 to 2015), and time of the day. For this
research electron densities were extracted from the IRI based on the following inputs:

Required input parameters
Year= 2000., Month= 03, Day= 01, Hour=1.5,
Time_type = Universal
Coordinate_type = Geographic
Latitude= 50., Longitude= 40., Height= 100.
Prof. parameters: Start= 60. Stop= 220. Step= 1.

Optional input parameters:
Sunspot number(Rz12) =not specifyed
Ionospheric index(IG12) =not specifyed
Upper limit for Electron content = not specifyed
F peak model = URSI
Ne Topside = NeQuick
foF2 Storm model = on
Bottomside Thickness = B0 Table
F1 occurrence probability = Scotto-1997 no L
D-Region Ne = IRI-95
Topside Te = TTSA-2000
Ion Composition = DS95/TT05
A value of -1 indicates that the parameter is not available for the specified range
TEC=-1, means you have not entered an upper boundary height in the OPTIONAL INPUT section.

Abstract (if available)

Abstract One of the natural sources of the electric current is the lightning discharge in the ionosphere. The electro-magnetic wave produced by intense lightning, like sprites, can significantly alter electrostatic, electromagnetic, acoustic, radar, and radio-frequency measurements in ionosphere and on the ground. There are numerous studies performed on the sprites and other types of lightning discharge activities. This study concentrates on very low frequency (VLF, 3-20 kHz) and extremely low frequency (ELF, 3-3000 Hz) electromagnetic waves radiated by lightning currents. These frequencies have very low attenuation rate and make the ground a good enough conductor where the radio wave arrives at the surface entering the ground nearly independent of the incident angle. The Earth curvature can be neglected for frequencies up to ~12 kHz since the wave incident angles are very small. Method being used is Finite-Difference Time-Domain (FDTD) in the two-dimensional cylindrical simulation environment. Surface impedance boundary condition is used for the air-ground interface

Linked assets

University of Southern California Dissertations and Theses