A Stochastic Finite-Difference Time-Domain (FDTD) Method for Assessing Material and Geometric Uncertainties in Rectangular Objects ^†

Abstract

The uncertainties present in a variety of electromagnetic (EM) problems may have important effects on the output parameters of interest. Unfortunately, deterministic schemes are not applicable in such cases, as they only utilize the nominal value of each random variable. In this work, a two-dimensional (2D) finite-difference time-domain (FDTD) algorithm is presented, which is suitable for assessing randomness in the electrical properties, as well as in the dimensions of orthogonal objects. The proposed technique is based on the stochastic FDTD method and manages to extract the mean and the standard deviation of the involved field quantities in one realization. This approach is applied to three test cases, where uncertainty exists in the electrical and geometrical parameters of various materials. The numerical results demonstrate the validity of our scheme, as similar outcomes are extracted compared to the Monte Carlo (MC) algorithm.

1. Introduction

Wave propagation in a material with random electric and geometric properties may introduce uncertainty in the involved field quantities. For example, the variability present in biologic tissues [1,2] may have a non-trivial effect on the absorbed radio wave power. Consequently, a reliable assessment of such randomness is an essential step, since neglecting the latter may increase the risk of misinterpretation of the simulated outcomes. Unfortunately, deterministic analysis is not suitable in those cases, as it only takes into account the nominal values of each random variable. For this reason, several methods exist in the literature for solving stochastic electromagnetic (EM) problems.

The most common technique for uncertainty assessment is the Monte Carlo (MC) algorithm [3]. Particularly, this approach performs a high number of sample realizations, each time selecting a different value of each random input. Then, the extracted data are post-processed for the computation of the statistical moments of the examined outputs. However, a major disadvantage of this scheme is the slow rate of convergence $O\left(1/\sqrt{N}\right)$, with N denoting the number of samples. As a result, this method requires elongated computational times for solving complex EM problems, thus rendering it impractical in such scenarios.

An alternative surrogate is the generalized polynomial chaos (PC) [4], which has been applied in a variety of uncertainty problems. Specifically, this technique is based on the orthogonal expansion of a second-order random process and can be more efficient than the MC algorithm, for cases with a low number of random variables. Unfortunately, the PC scheme suffers from the ‘curse of dimensionality‘, which means that the computational cost of this approach increases exponentially with the number of stochastic inputs or with the order of the expansion.

A more recent scheme combines the finite-difference time-domain (FDTD) algorithm with the Delta method [5], aiming at speeding up the stochastic analysis of random media problems. The proposed technique is named stochastic FDTD (S-FDTD) and manages to compute the mean and the standard deviation of the field components in one simulation. Therefore, the S-FDTD algorithm can be less computationally demanding than both the MC and the PC approaches. However, the accuracy of this scheme depends on correlation coefficient values, which appear in the updated equations.

In our work, we introduce an extension of the S-FDTD technique, which is able to assess random media, as well as uncertainties in the dimensions of rectangular dielectrics. The proposed technique enhances the work in [6] by additionally taking into consideration the uncertainties present in the conductivity of the random medium. The validation of our algorithm is performed via three two-dimensional (2D) stochastic EM problems. Numerical results prove the effectiveness of our approach, and comparisons with the MC method show that satisfactory outcomes can be extracted with a lower computational cost.

2. Literature Review of Related Works

The S-FDTD technique has been utilized in a variety of random media uncertainty scenarios. For example, the works in [1,5] introduced the mathematical background of this scheme and applied it in one dimensional (1D) biolectromagnetic problems. Moreover, the authors in [7] presented an extension of the S-FDTD algorithm, which is able to assess uncertainties caused by the electron and the ion content of the atmosphere. Their numerical outcomes indicated a significant efficiency enhancement compared to the MC technique, thus proving that the proposed method can be a good alternative for EM ionospheric propagation studies. A variation of the S-FDTD technique for assessing uncertainties in cold plasma medium is examined in [8]. The suggested approach is based on the bilinear transform technique and has been successfully tested in 1D, 2D, and three-dimensional (3D) cases.

The work in [9] introduced a novel S-FDTD scheme, which is free from time-step stability restrictions. This method utilizes the locally one-dimensional approach (hence named S-LOD-FDTD) and is useful for cases with dense grids. Numerical outcomes in 1D transmission-line examples indicate that reliable results can be extracted more efficiently than with the MC or the traditional S-FDTD algorithms. A further extension of the S-LOD FDTD technique to 2D EM uncertainty problems was presented in [10]. Larger computational time-steps render the proposed scheme faster than the MC and the S-FDTD approaches, in this case as well. The S-FDTD algorithm has also been combined with the alternating-direction implicit (ADI) method and applied to 1D cases with random dielectric materials [11]. Furthermore, in [12], a variation of the S-FDTD technique was presented, which is suitable for assessing randomness in circular dielectrics. Specifically, the Delta method was applied to the polar FDTD scheme and more reliable outcomes than the conventional S-FDTD approach are extracted for those problems.

In [13], a stochastic version of the finite-difference frequency-domain algorithm was proposed. This approach was tested on 1D and 2D bioelectromagnetic problems and manages to calculate the mean and the variance of the involved field quantities in one realization with a fair degree of accuracy. The work in [14] presents a modified S-FDTD technique, which is appropriate for assessing uncertainties in dispersive materials. Particularly, the proposed scheme is based on a D-H formulation and can be applied to general dispersive media of a single pole. Finally, the authors in [15] suggested an S-FDTD method that is able to assess geometrical variations in bioelectromagnetic problems. Their approach focused on the variance calculation of the specific absorption rate (SAR) in a 2D slice of a human model, as well as in estimating the impact of tissue variation in 1D and 2D biological cases.

3. Methodology

According to [16], the equations for the 2D FDTD algorithm in a lossy medium are:

$$\begin{array}{cc}\hfill & E_x{|}_{i+\frac{1}{2},j}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n=\frac{2\Delta t}{\left(2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t\right)\Delta y}\left(H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}-H_z{|}_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}\right)\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & E_y{|}_{i,j+\frac{1}{2}}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_y{|}_{i,j+\frac{1}{2}}^n=-\frac{2\Delta t}{\left(2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t\right)\Delta x}\left(H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}-H_z{|}_{i-\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right)\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}=\frac{\Delta t}{\mathsf{\mu }}\left(\frac{E_x{|}_{i+\frac{1}{2},j+1}^{n+1}-E_x{|}_{i+\frac{1}{2},j}^{n+1}}{\Delta y}-\right.\left.\frac{E_y{|}_{i+1,j+\frac{1}{2}}^{n+1}-E_y{|}_{i,j+\frac{1}{2}}^{n+1}}{\Delta x}\right)\hfill \end{array}$$

(3)

where $\mathsf{ϵ}$, $\mathsf{\mu }$, and $\underline{\mathsf{\sigma }}$ denote the dielectric permittivity, the magnetic permeability, and the conductivity of the background medium, respectively. Additionally, $\Delta x$ and $\Delta y$ are the spatial steps in the x and y directions, and $\Delta t$ is the time-step value. A medium with stochastic electric and magnetic properties is modeled by introducing uncertainty in the parameters $\mathsf{ϵ}$, $\mathsf{\mu }$, and $\underline{\mathsf{\sigma }}$. Also, the randomness in the width and height of the random dielectric is incorporated in $\Delta x$ and $\Delta y$. In our work, the magnetic permeability is considered deterministic, however this quantity can be also treated as stochastic by following a similar procedure to the one demonstrated below. As described in [5], the mean and the variance of a function f that depends on the variables $x_1$, $x_2$, ⋯, $x_n$ can be approximated by the following equations respectively:

$$\begin{array}{c}\hfill M\left\{f\left(x_1,x_2,\dots ,x_n\right)\right\}\approx f\left(M\left\{x_1\right\},M\left\{x_2\right\},\dots ,M\left\{x_n\right\}\right)\end{array}$$

(4)

$$\begin{array}{c}\hfill {\sigma }^2\left\{f(x_1,\dots ,x_n)\right\}\approx \sum _{i=1}^n\sum _{j=1}^n{\left.\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\right|}_{M\left\{x_1\right\},\dots ,M\left\{x_n\right\}}M\left\{(x_i-M\left\{x_i\right\})(x_j-M\left\{x_j\right\})\right\}\end{array}$$

(5)

where $M\left\{x_1\right\},M\left\{x_2\right\},\dots ,M\left\{x_n\right\}$ are the mean values of $x_1$, $x_2$, ⋯, $x_n$. According to (4), the mean field updates can be calculated directly from the standard 2D FDTD equations, as long as the random variables are substituted by their corresponding nominal values. In order to approximate the standard deviation of the field components, we start by applying the variance operator to both sides of (1):

$${\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}={\sigma }^2\left\{\frac{2\Delta t}{\left(2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t\right)\Delta y}\left(H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}-H_z{|}_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}\right)\right\}.$$

(6)

Next, we will continue with the expansion of the left side of (6). Utilizing the formula, ${\sigma }^2\left\{aX\pm bY\right\}=a^2{\sigma }^2\left\{X\right\}+b^2{\sigma }^2\left\{y\right\}\pm 2abCov\left\{X,Y\right\}$, where $Cov\left\{X,Y\right\}$ denotes the covariance between the random variables X and Y, yields:

$$\begin{array}{c}\hfill {\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n{E}_x{|}_{i+\frac{1}{2},j}^n\right\}={\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}+{\sigma }^2\left\{\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}\\ \hfill -2Cov\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1},\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}.\end{array}$$

(7)

Applying the identity $Cov\left\{X,Y\right\}={\rho }_{X,Y}\sigma \left\{X\right\}\sigma \left\{Y\right\}$, where ${\rho }_{X,Y}$ is the correlation coefficient between the variables X and Y, to (7) gives:

$$\begin{array}{c}\hfill {\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}={\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}+{\sigma }^2\left\{\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}\\ \hfill -2{\rho }_{E,\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}E}\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}\sigma \left\{\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}.\end{array}$$

(8)

The work in [5] suggests that ${\rho }_{E,\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}E}$ can be safely set equal to 1, thus leading to:

$${\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}={\left(\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}-\sigma \left\{\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}\right)}^2.$$

(9)

The electric field component is correlated with the conductivity and the dielectric permittivity. In this case, the approximation of (5) can be applied, as:

$$\begin{array}{c}\hfill {\sigma }^2\left\{\underset{g}{\underbrace{\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n}}\right\}\approx {\left(\frac{\partial g}{\partial E_x{|}_{i+\frac{1}{2},j}^n}\right)}^2{\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}+{\left(\frac{\partial g}{\partial \mathsf{ϵ}}\right)}^2{\sigma }^2\left\{\mathsf{ϵ}\right\}+{\left(\frac{\partial g}{\partial \underline{\mathsf{\sigma }}}\right)}^2{\sigma }^2\left\{\underline{\mathsf{\sigma }}\right\}+\cdots \\ \hfill +2\frac{\partial g}{\partial E_x{|}_{i+\frac{1}{2},j}^n}\frac{\partial g}{\partial \mathsf{ϵ}}\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}\sigma \left\{\mathsf{ϵ}\right\}{\rho }_{\mathsf{ϵ},E}+\cdots .\end{array}$$

(10)

Taking the square root of (10) and performing the approximations presented in [5], yields:

$$\begin{array}{c}\hfill \sigma \left\{\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}\approx \frac{2M\left\{\mathsf{ϵ}\right\}-\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}}{2M\left\{\mathsf{ϵ}\right\}+\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}}\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}\\ \hfill +\frac{4\Delta t\left(M\left\{\underline{\mathsf{\sigma }}\right\}{\rho }_{\mathsf{ϵ},E}\sigma \left\{\mathsf{ϵ}\right\}-M\left\{\mathsf{ϵ}\right\}{\rho }_{\underline{\mathsf{\sigma }},E}\sigma \left\{\underline{\mathsf{\sigma }}\right\}\right)}{{\left(2M\left\{\mathsf{ϵ}\right\}+\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}\right)}^2}M\left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}.\end{array}$$

(11)

As a result, (9) becomes:

$$\begin{array}{c}\hfill {\sigma }^2\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}-\frac{2\mathsf{ϵ}-\underline{\mathsf{\sigma }}\Delta t}{2\mathsf{ϵ}+\underline{\mathsf{\sigma }}\Delta t}{E}_x{|}_{i+\frac{1}{2},j}^n\right\}\approx \left(\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}-\frac{2M\left\{\mathsf{ϵ}\right\}-\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}}{2M\left\{\mathsf{ϵ}\right\}+\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}}\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}\right.\\ \hfill {\left.-\frac{4\Delta t\left(M\left\{\underline{\mathsf{\sigma }}\right\}{\rho }_{\mathsf{ϵ},E}\sigma \left\{\mathsf{ϵ}\right\}-M\left\{\mathsf{ϵ}\right\}{\rho }_{\underline{\mathsf{\sigma }},E}\sigma \left\{\underline{\mathsf{\sigma }}\right\}\right)}{{\left(2M\left\{\mathsf{ϵ}\right\}+\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}\right)}^2}M\left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}\right)}^2.\end{array}$$

(12)

The rest of the mathematical analysis is presented in Appendix A. The extracted equations are shown below:

$$\begin{array}{c}\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}\approx C_1\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}+C_{2}M\left\{E_x{|}_{i+\frac{1}{2},j}^n\right\}+C_3\left(\frac{\sigma \left\{H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}-\left\{H_z{|}_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}\right\}}{M\left\{\Delta y\right\}}\right)\end{array}$$

$$\begin{array}{c}-C_3\left(C_4+C_5\right)\left(M\left\{H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}-M\left\{H_z{|}_{i+\frac{1}{2},j-\frac{1}{2}}^{n+\frac{1}{2}}\right\}\right)\\ \sigma \left\{E_y{|}_{i,j+\frac{1}{2}}^{n+1}\right\}\approx C_1\sigma \left\{E_y{|}_{i,j+\frac{1}{2}}^{n+1}\right\}+C_{2}M\left\{E_y{|}_{i,j+\frac{1}{2}}^{n+1}\right\}-C_3\frac{\sigma \left\{H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}-\left\{H_z{|}_{i-\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}}{M\left\{\Delta x\right\}}\end{array}$$

(13)

$$\begin{array}{c}+C_3\left(C_4+C_6\right)\left(M\left\{H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}-M\left\{H_z{|}_{i-\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}\right)\end{array}$$

(14)

$$\begin{array}{c}\sigma \left\{H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}\right\}\approx \sigma \left\{H_z{|}_{i+\frac{1}{2},j+\frac{1}{2}}^{n-\frac{1}{2}}\right\}+D_1\left(\frac{\sigma \left\{E_x{|}_{i+\frac{1}{2},j+1}^{n+1}\right\}-\sigma \left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}}{M\left\{\Delta y\right\}}\right.\\ \left.-\frac{\sigma \left\{E_y{|}_{i+1,j+\frac{1}{2}}^{n+1}\right\}-\sigma \left\{E_y{|}_{i,j+\frac{1}{2}}^{n+1}\right\}}{M\left\{\Delta x\right\}}\right)+D_1{D}_2\left(M\left\{E_x{|}_{i+\frac{1}{2},j+1}^{n+1}\right\}-M\left\{E_x{|}_{i+\frac{1}{2},j}^{n+1}\right\}\right)\\ +D_1{D}_3\left(M\left\{E_y{|}_{i+1,j+\frac{1}{2}}^{n+1}\right\}-M\left\{E_y{|}_{i,j+\frac{1}{2}}^{n+1}\right\}\right)\end{array}$$

(15)

where

$$\begin{array}{c}{C}_1=\frac{2M\left\{\mathsf{ϵ}\right\}-\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}}{2M\left\{\mathsf{ϵ}\right\}+\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}},C_2=\frac{4\Delta t\left(M\left\{\underline{\mathsf{\sigma }}\right\}{\rho }_{\mathsf{ϵ},E}\sigma \left\{\mathsf{ϵ}\right\}-M\left\{\mathsf{ϵ}\right\}{\rho }_{\underline{\mathsf{\sigma }},E}\sigma \left\{\underline{\mathsf{\sigma }}\right\}\right)}{{\left(2M\left\{\mathsf{ϵ}\right\}+\Delta tM\left\{\underline{\mathsf{\sigma }}\right\}\right)}^2},\\ C_3=\frac{2\Delta t}{2M\left\{\mathsf{ϵ}\right\}+M\left\{\underline{\mathsf{\sigma }}\right\}\Delta t},C_4=\frac{2\sigma \left\{\mathsf{ϵ}\right\}{\rho }_{\mathsf{ϵ},H}+\Delta t\sigma \left\{\underline{\mathsf{\sigma }}\right\}{\rho }_{\underline{\mathsf{\sigma }},H}}{2M\left\{\mathsf{ϵ}\right\}+M\left\{\underline{\mathsf{\sigma }}\right\}\Delta t},C_5=\frac{\sigma \left\{\Delta y\right\}{\rho }_{\Delta y,H}}{M\left\{\Delta y\right\}}\\ C_6=\frac{\sigma \left\{\Delta x\right\}{\rho }_{\Delta x,H}}{M\left\{\Delta x\right\}}{D}_1=\frac{2\Delta t}{{\mathsf{\mu }}_0},D_2=\frac{{\rho }_{\Delta y,E}\sigma \left\{\Delta y\right\}}{M^2\left\{\Delta y\right\}},D_3=\frac{{\rho }_{\Delta x,E}\sigma \left\{\Delta x\right\}}{M^2\left\{\Delta x\right\}}.\end{array}$$

(16)

4. Numerical Results

The proposed technique is evaluated via three 2D EM problems with random dielectric slabs. Specifically, the first test case involves a lossless rectangular material with uncertainty in its width and height. In Figure 1, the schematic of the examined problem is presented. The computational domain is discretized into $200\times 200$ cells with mean spatial steps equal to ${\lambda }_0/60$, where ${\lambda }_0$ denotes the wavelength in the vacuum at 2 GHz. The random dielectric is located at the 140th cell, while the remaining space is filled with air. The dielectric permittivity of the slab has a value of $\mathsf{ϵ}=8{\mathsf{ϵ}}_0$ and the stochastic variables follow a normal distribution with standard deviations equal to the $5\%$ of their corresponding mean values. It is worth mentioning that the time-step is equal to $\Delta t=\frac{1}{2c_0\sqrt{\frac{1}{M^2\left\{\Delta x\right\}}+\frac{1}{M^2\left\{\Delta y\right\}}}}$ and the total simulation time is 7.2 ns. Finally, unwanted reflections are reduced by employing the first-order Mur absorbing boundary [17] condition (ABC). This ABC is simple and requires less computational resources than the perfectly matched layer (PML) [18]. The proposed S-FDTD technique is compared with 1000 MC realizations. Particularly, Figure 2 and Figure 3 depict the mean and the standard deviation of the magnetic field at point (140, 140). In Figure 3, the accuracy of the proposed technique is examined with various correlation coefficient values. By setting $\rho =0.8$, we guarantee a good agreement from 4 to 5 ns compared to the MC method. For $\rho =0.9$, valid outcomes can be extracted around 0 to 4 ns, while $\rho <0.8$ gives satisfactory results, between 5.25 and 7.2 ns.

The second scenario exhibits the same problem setup as in the first case, however an additional uncertainty in the dielectric permittivity of the medium is incorporated. Again, each random variable follows a normal distribution with standard deviation equal to the $5\%$ of the corresponding nominal value. In Figure 4 and Figure 5, the mean and the standard deviation of the $H_z$ component are depicted at point (120, 140). Accurate results compared to the MC technique are obtained by setting $\rho =0.5$.

The third numerical test exhibits two rectangular materials with randomness present in their dielectric permittivities, conductivities, and in their dimensions. Figure 6 illustrates the schematic of this scenario, where the mean dielectric permittivity of the top material is equal to $M\left\{\mathsf{ϵ}\right\}=8{\mathsf{ϵ}}_0$ and the corresponding conductivity has a value $M\left\{\underline{\mathsf{\sigma }}\right\}=0.2$ S/m. Similarly, the same parameters for the bottom dielectric are $M\left\{\mathsf{ϵ}\right\}=6.2{ϵ}_0$ and $M\left\{\underline{\mathsf{\sigma }}\right\}=0.4$ S/m. Also, the variability of each stochastic input is the same as in the case of the first two scenarios. Figure 7 depicts the mean $H_z$ field component, while Figure 8 illustrates the standard deviation of the same quantity. In this test case, satisfactory outcomes are extracted by setting $\rho$ to be around 0.7. It is also worth mentioning that for all problems, the required simulation time for the MC technique is approximately 1.44 h, while the proposed S-FDTD approach needs 17.88 s. Unfortunately, the choice of correlation coefficient values is still an open problem, as they cannot be known beforehand. Setting $\rho =1$ can often to lead to an overestimation of the standard deviation, which can be observed in Figure 3, Figure 5 and Figure 8. According to [19], a partial solution to this problem is to perform a small fraction of MC simulations, which unavoidably increases the total computational cost of the S-FDTD technique.

5. Conclusions

An extension of the S-FDTD algorithm is presented, which is able to assess uncertainties in the electrical and geometrical characteristics of rectangular dielectrics. The proposed approach manages to extract the mean and the standard deviation of the involved field quantities in one realization, thus rendering it less computationally expensive than the MC method. Nevertheless, the accuracy of this scheme highly depends on the correlation coefficient values, which unfortunately cannot be known beforehand. Particularly, those coefficients are coupled with the random field components and, therefore, their estimation is infeasible without the statistical characteristics of the latter. In future works, the proposed approach will be extended to 3D cases, by applying the stochastic analysis to the 3D FDTD equations. As a result, geometric uncertainties presented in the length, width, and height of random slabs can be assessed. In this case the S-FDTD approach will consist of twelve updated equations, as the mean and the standard deviation of the six field components need to be computed in those scenarios.