Wave Equation FDTD
It is possible to use the FDTD method using the second-order wave equation version of Maxwell equations, instead of the usual first order curl equations (The standard Yee algorithm).
FDTD using the wave equation was first described in Aoyagi1993. This article builds upon the previous work of Aoyagi etal (Aoyagi1993) and includes a “wave equation” version of Berenger’s PML ABC.
It is also adapted from the earlier standard 2-D Curl Equation code of Dr. Susan C. Hagness (See Taflove2000)
The FDTD wave equation formulation probably has no practical value, as implementing the wave equation version of the PML is much more memory intensive and compute intensive than the standard Yee algorithm.
The point of the wave equation is to “decouple” the H and E fields, but the penalty is that it is then necessary to keep around an old value of Ex and Ey for every point on the grid.
This program is essentially a “novelty” and is provided as an existence proof. Implementing the PML as a wave equation can be done. (But you probably wouldn’t want to)
The Wave Equation form of the Maxwell’s Equations
Starting from the modified Maxwell’s equations for FDTD (Equation 2.10 and Equation 2.9):
Taking the partial time derivative of Equation 4.1 gives:
And substituting Equation 4.2 into Equation 4.3 gives:
Re-Arranging Equation 4.4 gives:
where,
 |
It is difficult to cleanly reduce the equation much further unless the magnetic resistivity is a constant (for all practical purposes, 0.)
For materials such as Berenger’s PML, the magnetic resistivity will not be a constant.
However, for cases when the magnetic resistivity is 0, then Equation 4.5 reduces to:
Now given the vector identity Equation C.20:
Substitute Equation C.20 (Equation 4.7) into Equation 4.6 gives:
Wave Equation (with losses and sources)
Which is the “Full Wave Equation” for E, including sources and losses.
Wave Equation (for free space)
In a lossless region with no sources Equation 4.8 reduces to the standard wave equation for free space:
When…
 |
Then,
The 2-D TEz Wave Equation FDTD Algorithm
The FDTD Wave Equation algorithm
solves for electric and magnetic fields using the uncoupled (E, H are uncoupled)
second order wave equations. (As opposed to solving for both the electric or magnetic field using the first order coupled Maxwell’s curl equations)
For the 2-D wave equation the components of E (Ex,Ey) are uncoupled from the components of H (Hz).
The wave equation algorithm uses the same Yee grid as the standard curl equation algorithm:
Close-up view of the Yee Grid for the 2-D TEz (transverse electric mode)
Generating the FDTD Wave Equation algorithm
It may be easier to derive the wave equation by starting with the
difference equations for the standard curl equations algorithm:
where,
caex[i][j] = (1.0 - temporaryX) / (1.0 + temporaryX) cbex[i][j] = (dt / (permittivityX[i][j] * dx)) / (1.0 + temporaryX) caey[i][j] = (1.0 - temporaryY) / (1.0 + temporaryY) cbey[i][j] = (dt / (permittivityY[i][j] * dx)) / (1.0 + temporaryY) dahz[i][j] = (1.0 - temporaryZ) / (1.0 + temporaryZ) dbhz[i][j] = (dt / (permeability[i][j] * dx) / (1.0 + temporaryZ) and, temporaryX = dt * conductivityX[i][j] / (2.0 * permittivityX[i][j] ) temporaryY = dt * conductivityY[i][j] / (2.0 * permittivityY[i][j] ) temporaryZ = dt * magneticResistivity[i][j] / (2.0 * permeability[i][j]) from Taflove1995 p.67
Example: Generating the FDTD Wave Equation for Ex
Starting from Equation 4.11:
the equation for Ex at time ‘n’ is:
Subtracting Equation 4.14 from Equation 4.13 gives:
Now assume that magneticResistivity[i][j] = 0 and permeability[i][j] is a constant, then dahz and dbhz reduce to:
dahz = 1.0 dbhz = (dt / (permeability * dx)
Given that dahz = 1.0 and dbhz is a constant, Equation 4.10 can be rewritten as:
Likewise,
Now, substituting Equation 4.16 and Equation 4.17 into Equation 4.15 gives:
Which can be reduced to:
Which is the wave equation for Ex in a region with no sources and
magneticResistivity[i][j] = 0 and permeability[i][j] is a constant.
A similar procedure can be followed to obtain Ey.
FDTD Wave Equation in a region with a Source
In the following equations, i=xsource, j=ysource…
The following implements a Hz Magnetic infinite line source at x=i+.5 and y=j+.5
To implement the wave equation formulation without using Hz, equations for the four Ex, Ey, surrounding
Hz(i+.5,j+.5) must be created. These are: Ex(i,j), Ex(i,j+1), Ey(i,j), and Ey(i+1,j).
Example: Generating the Source FDTD Wave Equation for Ex at x=i, y=j
Starting from Equation 4.15:
Now assume that conductivity[i][j] = 0 and permittivity[i][j] is a constant, magneticResistivity[i][j] = 0 and permeability[i][j] is a constant, then caex, cbex, caey, cbey, dahz and dbhz reduce to:
caex = caey = 1.0 cbex = cbey = (dt / (permittivity * dx) dahz = 1.0 dbhz = (dt / (permeability * dx)
Therefore, Equation 4.20 reduces to:
Now, substituting Equation 4.17 into Equation 4.21 gives:
Which is the source equation for Ex(i,j).
By a similar procedure, Ex(i,j+1), Ey(i,j), Ey(i+1,j) can be obtained.
The actual equations are shown below:
i = xSource ;
j = ySource;
temporary = cbex * (currentSourceValue - oldSourceValue);
exNew[i][j] = (2.0 - cbex * dbhz) * ex[i][j] - exOld[i][j] + temporary
+ (cbex * dbhz) * ( ex[i][j-1] - ey[i][j-1] + ey[i+1][j-1] )
exNew[i][j+1] = (2.0 - cbex * dbhz) * ex[i][j+1] - exOld[i][j+1] - temporary
+ (cbex * dbhz) * ( ex[i][j+2] + ey[i][j+1] - ey[i+1][j+1] )
eyNew[i][j] = (2.0 - cbex * dbhz) * ey[i][j] - eyOld[i][j] - temporary
+ (cbex * dbhz) * ( ey[i-1][j] - ex[i-1][j] + ex[i-1][j+1] )
eyNew[i+1][j] = (2.0 - cbex * dbhz) * ey[i+1][j] - eyOld[i+1][j] + temporary
+ (cbex * dbhz) * ( ey[i+2][j] + ex[i+1][j] - ex[i+1][j+1] )
These equations simulate a Magnetic infinite line source Hz(i+.5,j+.5). Note: Trying to simulate an Hz hertzian dipole using the wave equation implementation is much more complicated compared to the standard Yee curl equations implementation!
Berenger Perfectly Matched Layer (PML) Absorbing Boundary Conditions (ABCs)
The Berenger PML can be implemented using wave equations.
Starting with the split-field difference equations for the standard curl equations algorithm:
where,
caex[i][j] = (1.0 - temporaryX) / (1.0 + temporaryX) cbex[i][j] = (dt / (permittivity0 * dx)) / (1.0 + temporaryX) caey[i][j] = (1.0 - temporaryY) / (1.0 + temporaryY) cbey[i][j] = (dt / (permittivity0 * dx)) / (1.0 + temporaryY) dahzx[i+.5][j+.5] = (1.0 - temporaryZx) / (1.0 + temporaryZx) dbhzx[i+.5][j+.5] = (dt / (permeability0 * dx) / (1.0 + temporaryZx) dahzy[i+.5][j+.5] = (1.0 - temporaryZy) / (1.0 + temporaryZy) dbhzy[i+.5][j+.5] = (dt / (permeability0 * dx) / (1.0 + temporaryZy) and, temporaryX = dt * conductivityY[i][j] / (2.0 * permittivity0 ) temporaryY = dt * conductivityX[i][j] / (2.0 * permittivity0 ) temporaryZx = dt * magneticResistivityX[i+.5][j+.5] / (2.0 * permeability0) temporaryZy = dt * magneticResistivityY[i+.5][j+.5] / (2.0 * permeability0) from Taflove1995 p.183
In the 2-D TEz FDTD split-field PML curl equations, only Hz is split (Hzx, Hzy). In the wave equation version, Hz is not used,
so Ex and Ey must be split. It was found that the following Ex and Ey split-fields correctly modeled the Berenger PML:
Where,
And
Example: Generating the PML FDTD Wave Equation for Exxa
Starting from Equation 4.27, the equation for exxa(i,j,n) is:
Now multiply Equation 4.37 by dahzx(i+.5,j+.5), to get:
Now subtract Equation 4.38 from Equation 4.27, to get:
Now, substituting Equation 4.25 into Equation 4.39 and rearranging gives:
Where Equation 4.40 is the split-field PML equation for exxa.
By a similar procedure, equations for exxb, exya, exyb, eyxa, eyxb, eyya, eyyb can be obtained.
The actual equations are shown below:
exxaNew[i][j] = (caex[i][j] + dahzx[i][j]) * exxa[i][j]
- (caex[i][j] * dahzx[i][j]) * exxaOld[i][j]
+ (cbex[i][j] * dbhzx[i][j]) * (ey[i][j] - ey[i+1][j]);
exxbNew[i][j] = (caex[i][j] + dahzx[i][j-1]) * exxb[i][j]
- (caex[i][j] * dahzx[i][j-1]) * exxbOld[i][j]
- (cbex[i][j] * dbhzx[i][j-1]) * (ey[i][j-1] - ey[i+1][j-1]);
exyaNew[i][j] = (caex[i][j] + dahzy[i][j]) * exya[i][j]
- (caex[i][j] * dahzy[i][j]) * exyaOld[i][j]
+ (cbex[i][j] * dbhzy[i][j]) * (ex[i][j+1] - ex[i][j]);
exybNew[i][j] = (caex[i][j] + dahzy[i][j-1]) * exyb[i][j]
- (caex[i][j] * dahzy[i][j-1]) * exybOld[i][j]
- (cbex[i][j] * dbhzy[i][j-1]) * (ex[i][j] - ex[i][j-1]);
eyxaNew[i][j] = (caey[i][j] + dahzx[i-1][j]) * eyxa[i][j]
- (caey[i][j] * dahzx[i-1][j]) * eyxaOld[i][j]
+ (cbey[i][j] * dbhzx[i-1][j]) * (ey[i-1][j] - ey[i][j]);
eyxbNew[i][j] = (caey[i][j] + dahzx[i][j]) * eyxb[i][j]
- (caey[i][j] * dahzx[i][j]) * eyxbOld[i][j]
- (cbey[i][j] * dbhzx[i][j]) * (ey[i][j] - ey[i+1][j]);
eyyaNew[i][j] = (caey[i][j] + dahzy[i-1][j]) * eyya[i][j]
- (caey[i][j] * dahzy[i-1][j]) * eyyaOld[i][j]
+ (cbey[i][j] * dbhzy[i-1][j]) * (ex[i-1][j+1] - ex[i-1][j]);
eyybNew[i][j] = (caey[i][j] + dahzy[i][j]) * eyyb[i][j]
- (caey[i][j] * dahzy[i][j]) * eyybOld[i][j]
- (cbey[i][j] * dbhzy[i][j]) * (ex[i][j+1] - ex[i][j]);
As can be seen, the wave equation version of the PML is much more memory and compute intensive compared to the standard Berenger PML curl equations.
It should be noted that the coefficients caex, cbex, caey, cbey, dahzx, dbhzx, dahzy, dbhzy are exactly the
same as in the standard Berenger PML curl equations. The wave equation uses the same PML gradient as the curl equations.
So the gradient code can be shared between programs.
FDTD Source Code and Results
- Wave2d.c .c file (FDTD code)
- Results of simulation (Display Hz vs time) Note: Displays Hz in main grid only
See Also
- Fdtd Main Index Page (with example code)
- Maxwell’s Equations
- The Yee Algorithm
- 2D UPML Algorithm and Theory
- The 2D FDTD Wave Hz Algorithm
- Analytical Solution: Infinite Magnetic Line-source
- Numerical Methods
- Spherical and Cartesian coordinates
- Vector Algebra
- Fundamental and Derived Units
- Glossary
- FDTD References
- FDTD Resources
References
- Aoyagi, Lee, Mittra; “A Hybrid Yee Algorithm/Scalar-Wave Equation Approach”; IEEE Transactions on Microwave Theory and Techniques, vol. 41, no. 9, September 1993
External Links
- Aoyagi etal 1993 at IEEE.org ($$)