Abstract
The paper reviews the application of deterministic-stochastic models in some areas of computational electromagnetics. Namely, in certain problems there is an uncertainty in the input data set as some properties of a system are partly or entirely unknown. Thus, a simple stochastic collocation (SC) method is used to determine relevant statistics about given responses. The SC approach also provides the assessment of related confidence intervals in the set of calculated numerical results. The expansion of statistical output in terms of mean and variance over a polynomial basis, via SC method, is shown to be robust and efficient approach providing a satisfactory convergence rate. This review paper provides certain computational examples from the previous work by the authors illustrating successful application of SC technique in the areas of ground penetrating radar (GPR), human exposure to electromagnetic fields, and buried lines and grounding systems.
1. Introduction
Some areas in computational electromagnetics (CEM) suffer from uncertainties of the input parameters resulting in the uncertainties in the assessment of the related response. These problems could be overcome, to a certain extent, by an efficient combination of well-established deterministic electromagnetic models with certain stochastic methods.
Traditional methods rely upon statistical approaches such as brute Monte Carlo (MC) simulations and various sampling techniques like stratified sampling and Latin hypercube sampling (LHS). These methods are easy to implement and do not suffer from “curse of dimensionality” which means that the sample size does not depend on number of random variables. On the other hand the sample size needs to be very high (>100.000), and thus they exhibit very slow convergence.
Contrary to statistical approaches, the nonstatistical based techniques aim to represent the unknown stochastic solution as a function of random input variables. Among the various methods available in the literature, the spectral discretization based technique—generalized polynomial chaos (gPCE)—emerged as the most used approach in the stochastic CEM. The gPCE framework comprises stochastic Galerkin method (SGM) and stochastic collocation method (SCM) for solving stochastic equations [1]. The SGMs have been successfully used in recent years in the area of circuit uncertainty modeling both for Signal Integrity (SI) [2] and microwave applications [3] as well as in stochastic dosimetry [4, 5]. The intrusive nature of SGM implies more demanding implementation since a development of new codes is required. On the contrary, the nonintrusive nature of SCM enables the use of reliable deterministic models as black boxes in stochastic computations. Both approaches exhibit fast convergence and high accuracy under different conditions and a detailed comparison of their use in EMC simulation can be found in [6].
The combination of the nonintrusive, sampling based nature of Monte Carlo simulations with the polynomial approximation of output value, which is the characteristic of gPCE methods, made stochastic collocation one of the most researched and applied stochastic approaches [7–9]. The examples of successful coupling of SCM and its variations with different deterministic codes have been reported in stochastic dosimetry [10, 11], in the analysis of complex power distribution networks (PDN) [12], in the design of integrated circuits (ICs), microelectromechanical systems, (MEMSs), and photonic circuits [13, 14], in the simulation of EM-circuit systems [15, 16], and in the area of antenna modeling [17–19] and electromagnetic compatibility (EMC) of space applications [20]. This paper reviews the previous work of the authors pertaining to the use of SC techniques in areas of CEM such as ground penetrating radar (GPR), EM-thermal dosimetry of the eye and brain, and buried lines and grounding systems [21]. Some illustrative computational examples for the transient transmitted field from GPR antenna [22, 23], specific absorption rate (SAR) distribution in the brain and the eye [24, 25], plane wave coupling to buried conductors [26, 27], and transient analysis of grounding electrodes [28] pertaining to certain statistical moments are given in the paper as well.
The paper is organized, as follows: first an outline of the SC method is given in Section 2 along with the short overview of its applications in various areas. Section 3 deals with different applications of SCM carried out by authors with related examples. This is followed by some conclusions and guidelines for a future work.
2. An Outline of the Stochastic Collocation Method
The uncertainty quantification (UQ) of the unknown stochastic output of the model is preceded by two steps: the UQ of input parameters and uncertainty propagation (UP) of uncertainties present in the model inputs to the output of interest. The UQ of input parameters implies modeling the input parameters as random variables and/or random processes. The UP refers to the choice and implementation of the stochastic method that is capable of solving the stochastic model. The advantage of the stochastic collocation method (SCM) used in this work is its simplicity, a strong mathematical background, and the polynomial representation of stochastic output. The nonintrusive nature of the method enables the use of deterministic models as a black box. This way, previously validated computational models, such as FEM-BEM models described in Section 3, are used at predetermined set of simulation points. This section outlines the fundaments of the mathematical background for the SCM, with the brief mention of some other variants.
Once the deterministic modeling of a problem of interest is completed, a stochastic processing of the numerical results can be carried out via the SC method [1, 22–27]. The theoretical basis of SC technique is to use the polynomial approximation of the considered output for a certain number of random parameters.
Without loss of generality, the theoretical principles can be presented taking into account a single random parameter . Within this framework, a given mapping through function can be assumed and a random variable (RV) given by an a priori probabilistic distribution ( standing for the relation between the output and input parameter, resp.). Consequently, the expressions for statistical moments, such as mean and variance of RV, can be derived.
First, an approximation of function over th order Lagrange polynomial basis can be written as follows:where Lagrange polynomials are given by
And one has the following property:where denotes Kronecker symbol.
It is worth mentioning that, although Lagrange polynomials are mostly used for the polynomial representation of the stochastic output, other types of basis functions are also possible. In [15] SCM is used to assess the uncertainty quantification for the hybrid EM-circuit systems. Two types of basis functions are used: Lagrange polynomials that have the character of locally global basis functions and piecewise multilinear basis functions that are able to capture the discontinuous issues in stochastic solutions.
Note that collocation points from (1) correspond to the points in Gauss quadrature rule attached to the probability distribution of random inputs (e.g., Legendre polynomials for a uniform distribution and Hermite polynomials for a Gaussian distribution). The choice of an adapted quadrature rule [19] yieldswhere is the probability density of random variable , while represents a function with a sufficient regularity and the coefficients are collocation weights chosen in a way to ensure an exact quadrature rule for polynomials with a degree lower or equal to .
However, the choice of collocation points is not limited to the points that supervene from Gauss quadrature rules. Other formulas may be used to generate the abscissas for the interpolation in (1) such as Clenshaw-Curtis formula with nonequidistant abscissas given as zeros of the extreme points of Chebyshev polynomials as in [15] or equidistant points in [14, 15]. These sets of points exhibit nested fashion unlike Gauss points, which is desirable in certain applications.
Higher dimensions yield multivariate interpolants in (1). The simplest way to interpolate is by using the tensor product in each random dimension. This approach is justified for the cases when random dimension is up to 5 [1]. Successful applications of a fully tensorized SCM model can be found in [7–9]. However, for random dimension > 5, SCM exhibits the “curse of dimensionality” and different techniques need to be considered. One of the most popular approaches is sparse grid interpolation based on Smolyak algorithm [1]. Sparse grid stochastic collocation method (SGSCM) has become an interest of many researchers. In [14] a conventional SGSCM has been applied for stochastic characterization of MEMS. An adaptive approach of SGSCM has been employed in [15] for the stochastic modeling of hybrid EM-circuits while a dimension-reduced sparse grid strategy has been proposed for SCM in EMC software in [16].
Furthermore, different cubature formulas for dealing with multidimensional integrals may be used for the generation of collocation points. One of the examples is Stroud cubature based SC method whose application for the variability analysis of PDNs within SPICE environments is reported in [12] and, for efficient computation of RCS from scatterers of uncertain shapes, in [19].
Once the desired output has been expressed as polynomial dependent on random input variables, the expression for the stochastic moments can be easily derived following their definitions from statistics. The expectation (mean) of is by definitionwhere is a given RV.
Next step is to replace in terms of Lagrange polynomials and expectation of can be written as follows:The Gauss quadrature rule yields
Taking into account the property of Lagrange polynomials (3) one obtains
Combining (6)–(8) leads to an expression for mean of :
Furthermore, the same approach is used to evaluate variance of , which is by definition given byIt can be written as
Again, utilizing the expansion of function over Lagrange polynomial basis (1) yieldswhich can be also written as
Furthermore, integral from the right-hand side of (13) according to Gauss quadrature rule (7) becomesFinally, inserting (14) into (13), one obtains and due to Kronecker property nonzero terms are only with leading to
Now the statistical moments around a considered output can be written as , where stands for a given order of statistical moment; for example, the mean of can then be written as
General expression is given as follows:where represents a mapping pertaining to a certain statistical moment , that is, with such as
Assuming the mutual independence of input random variables the generalization to multi-RVs case is straightforward. Instead of single random variable, one deals with a random vector and th statistical moment of output can be written as [28]where stands for weight of the expansion for random vector of size . Some more details related to this approach are available elsewhere, for example, in [28].
Although the assumption on independence of random input variables is valid in many CEM applications, it is worth noting that correlated random variables can also be taken into account. For example, in [29] collocation points that belong to correlated inputs are computed by using the different orthogonal polynomials, while in [30] the use of partitioning and clustering methods for generation of SC points is proposed which naturally extends to incorporation of correlated inputs.
3. Applications in Computational Electromagnetics
This section outlines the governing equations pertaining to CEM deterministic models reviewed in this work. Thus, first, Section 3.1 deals with transient field generated by GPR dipole antenna. Modeling of induced SAR in the brain and eye is presented in Section 3.2. Sections 3.3 and 3.4 deal with plane wave coupling to buried wire and current injection of horizontal grounding electrode, respectively.
3.1. Transient Field in the Soil Generated by GPR Dipole Antenna
The dipole is characterized by a length and radius and is located at height above an interface as depicted in Figure 1.
Direct time domain deterministic model of GPR dipole antenna is based on the space-time Hallen integral equation and corresponding integral formula for the calculation of the transmitted field into the lossy ground.
The transient current along the dipole antenna excited by an equivalent voltage generator is obtained by solving the Hallen integral equation [31]:where is velocity of light, is the free-space wave impedance, and and are the distances between observation point and source point on real and image wire, respectively. Time signals and account for the multiple reflections from the wire ends and can be determined by assuming the zero current at the wire ends [31]. The interface effects are taken into account via the reflection coefficient (RC) [31]:where
is the Dirac impulse; and is the th order modified Bessel function of the first kind. The Hallen equation (21) is solved via space-time version of Galerkin-Bubnov indirect boundary element method (GB-IBEM) [31].
The transient field transmitted into a lossy ground is given by an integral formula [31]where is velocity of wave propagation in the lower medium and is the distance from the dipole antenna to the observation point located in the lower medium:
The influence of the two-media interface is taken into account via the simplified space-time transmission coefficient arising from the modified image theory (MIT) [31]:where
The integral formula is handled via boundary element formalism as well. The mathematical details of these procedures could be found elsewhere, for example, in [31].
All computational examples are related to the horizontal dipole antenna of length = 1 m and radius = 6.74 mm above a real ground.
First, a simplified case of a dielectric half-space is considered [22]. The antenna is excited at its center by the Gaussian pulse voltage source:with parameters = 1 V, = s−1, = 1.43 ns.
Figure 2 shows the mean of the electric field transmitted in the soil for the following random variables: (antenna height above a soil), (soil permittivity), and (field observation point). First, one random variable at a time is considered, as follows:(i): central location, m; uniformly distributed: (in meters);(ii): initial value, ; uniformly distributed: ;(iii): initial value, = 1 m; uniformly distributed: (in meters).
Significant discrepancies appear between different scenarios with one random variable at a time, that is, between , , and models. The influence of each RV on the output is computed using the following relation:
Figure 3 shows the influence of each of three different RVs on the output .
Ranking random parameters from the most influential to the least influential (from -field variances), the model is more sensitive to RV3 (depth accuracy) and RV1 (antenna height) and less sensitive to soil permittivity.
Next example deals with the same GPR antenna radiating over a lossy ground with ground conductivity as the fourth random input variable (RV) [23]. The input RVs follow the uniform distribution. The observed types of soil are wet and average sand, respectively. The difference is in the corresponding ranges for relative permittivity of sand: and for wet and average sand, respectively. The other RVs are modeled as mS/m; cm; m.
Stochastic simulations are carried out for 4 univariate cases each with 3, 5, 7, and 9 sigma points and for 4-dimensional case with 54 input points. The impact factor is calculated according to (29).
Figure 4 depicts the mean value for the electric field transmitted into the soil versus time for wet sand when only permittivity is random. The relative permittivity is directly related to the propagation velocity; therefore, it has the highest impact on the time delay of a signal. More sigma points imply adding signals with different starting point to calculation. Therefore, for the given range of permittivity stochastic moments for the electric field versus time cannot be obtained. For this reason additional calculations were done for the case of average type of sand with narrower permittivity range.
Figure 5 shows the rank of input parameter for wet and average type of soil, respectively, with respect to transmitted electric field versus time. The relative permittivity has the highest impact followed by the depth of the observation point. The soil conductivity exhibits the minimal influence.
(a)
(b)
The maximum value of the transmitted field is mostly affected by the soil conductivity. The permittivity, the height, and the depth have a very similar influence on maximum field value for the wet sand. However, if the range of relative permittivity is reduced, the influence of relative permittivity is much smaller than the influence of the other two variables. Nevertheless, the multivariate simulation for both soil types has similar average value and standard deviation since the major origin of variations is the soil conductivity.
3.2. Induced SAR in the Brain and the Eye and Related Thermal Response
Deterministic frequency domain model of the human brain exposed to microwave radiation is based on the set of coupled electric field integral equations (EFIEs) and the related solution via method of moments (MoM) scheme. The mathematical details could be found elsewhere, for example, in [24, 32].
The main task of high frequency (HF) dosimetry is to quantify thermal effects, provided the distribution of the electromagnetic energy absorbed by the body is known. The power dissipated in the tissue is expressed in terms of the specific absorption rate (SAR) defined by the rate of energy absorbed by the unit body mass : where is the dissipated power, is the peak value of the electric field, respectively, is the tissue density, and σ is the tissue conductivity.
The temperature increase because the SAR induced in a tissue is obtained from the solution of the steady-state Pennes bioheat equation [24, 32]:where stands for the tissue thermal conductivity, is the blood mass density, is the specific heat capacity of blood, is the blood perfusion rate, is the arterial temperature, is the heat source due to metabolic processes, and term represents the volumetric heat source due to an external electromagnetic field.
The corresponding boundary condition at the interface between the model surface and the air iswhere is the heat flux density:while , , and are the convection coefficient, the temperature of the surface, and the temperature of the air, respectively.
The temperature rise in a tissue exposed to external HF fields is obtained by solving the bioheat transfer equation. The Pennes equation is solved via finite element method (FEM) and the details are available elsewhere, for example, in [24, 32].
The brain exposed to HF radiation is analyzed by solving the coupled surface integral equations (SIE) [32]:where and are the unknown equivalent electric and magnetic current densities, respectively, is the wave number of a medium , and is the interior/exterior Green function for the homogeneous medium [32] given bywhere is the distance from the source to the observation point, respectively, while denotes the index of an interior/exterior domain; that is, .
The set of integral equations (34) is solved via MoM procedure presented in [32].
Plane wave incident on the corneal part of the eye, treated as an exterior unbounded scattering problem, is formulated via the Stratton-Chu formulation; that is, the time-harmonic electric field at the exterior domain is expressed by the following boundary integral equation [25]:where is the incident electric field.
The interior domain is governed by the vector Helmholtz equation [25]:where subscripts and stand for the exterior and interior region, respectively.
This coupled formulation (36)-(37) is handled via the hybrid BEM/FEM technique.
The response of interest is the SAR distribution in the brain at = 900 MHz and = 1800 MHz, respectively, for the case of vertical polarization. The power density of the incident plane wave is mW/cm2.
The brain model is generated from a Google SketchUp rendering of the human brain. The dimensions of the average adult human brain are as follows: width 131.8 mm, length 161.1 mm, and height 139 mm, while the frequency dependent parameters of the human brain are presented in Table 1.
Figure 6 shows the results for the SAR and the resulting temperature increase obtained in the brain model at 900 MHz [24].
Each , where RVs are brain width, length, height, relative dielectric permittivity, and conductivity, respectively, is assumed to be uniformly distributed around deterministic values with common coefficient of variation equal to 5.77% [24].
Figure 7 shows the SC convergence with 3, 5, and 7 points for the SAR values calculations and different random variables ).
SC offers a precise assessment of the first statistical moment of maximum SAR with 5 multiphysics simulations, as the maximum value is between 0.87 and 0.9 W/kg, indicating the importance of modeling RVs at the same time.
Figures 8 and 9 are related to the results obtained from stochastic dosimetry modeling of the eye.
(a)
(b)
(a)
(b)
(c)
(d)
The deterministic eye parameter values are available elsewhere, for example, in [25], while only the relative permittivity and the electrical conductivity of vitreous body and cornea are considered here. The mean deterministic values and the stochastic variation assuming the uniform distribution, as the realistic distributions are not available, are given in Table 2.
Figure 8 shows the high convergence of the SC technique jointly with its huge efficiency as mostly only 3 or 5 simulations are required to compute mean and variance of the induced electric field without additional computational cost.
Figure 9 gives information of the first order sensitivity of the model to a corresponding random variable for the SAR assessment in the eye exposed to plane wave, that is, the SAR variance.
As evident from Figures 9(a) and 9(b), higher levels of SAR variance are obtained inside the central part of vitreous body while in the remaining parts this value is rather low. On the other hand, Figures 9(c) and 9(d) emphasizes the importance of the corneal part, as the distribution of SAR variance is mostly concentrated around the corneal region although some noticeable effects exist also in the anterior chamber.
3.3. Electromagnetic Field Coupling to Buried Conductors
The deterministic analysis of the plane wave scattering from buried wires is based on the space-time Pocklington integral equation. The geometry of a thin wire of length and radius is buried in a lossy medium at depth , excited by the plane wave as shown in Figure 10.
A transient response of a horizontal buried wire is governed by the following equation of the Pocklington type [26]:where is the unknown space-time dependent current, is the tangential component of the transmitted field, and is the space-time corresponding reflection coefficient arising from the modified image theory (MIT) approach, reported in [33]. Detailed derivation of (38) can be found in [34].
The distance from the source point in the wire axis to the observation point located on the wire surface iswhile the distance from the source point on the image wire to the observation point on the original wire, according to the image theory, is
Time constant and propagation velocity in the lossy medium are given with
The influence of the earth-air interface is taken into account via the reflection coefficient arising from the MIT and is given with [34]where the corresponding time constants are
Note that the reflection coefficient (42) represents rather simple characterization of the earth-air interface, taking into account only medium properties. An accuracy of (42) has been discussed in [34]. The solution is carried out analytically, as presented in [34].
The response of interest is the transient current at the center of the wire excited by a transmitted plane wave. The entire stochastic modeling is based upon realistic values of environmental parameters [26]:(i)soil conductivity : mS/m; uniformly distributed: (in mS/m);(ii)length of wire: m; uniformly distributed: (in meters);(iii)burial depth : m; uniformly distributed: (in meters).
Without loss of generality, the problem can be addressed following different assumptions about the statistical distribution laws.
Figure 11 shows mean (+one standard deviation) of the current at the wire center under uncertain constraints fully tensorized (i.e., with 3 RVs). The sensitivity analysis provides relevant information needed to decrease the total number of SC points required for each RV and optimize the “full-tensor” random model to an “asymmetrical” one.
Figure 12 provides convergence rates from the current variance including a complete random model: only 5 points are necessary to precisely describe the influence of random burying depth (RV3). Nearly zero levels of the current (mean and variance) below 0.03 μs involve instability of the convergence criterion (and positive SC gaps).
3.4. Transient Analysis of Grounding Electrode
This subsection deals with a transient analysis of a horizontal grounding electrode. The deterministic formulation is based on the homogeneous space-time Pocklington integral equation. A horizontal grounding electrode excited at one end with an equivalent current source of length and radius , buried in a lossy medium of conductivity and relative permittivity , at depth , is shown in Figure 13.
Thus, the current distribution along the wire is governed by the space-time homogeneous Pocklington integrodifferential equation [35]:The parameters used in (44) are already defined in Section 3.3.
Integrodifferential equation (44) is solved analytically, using certain approximations [35].
Additionally, transient impedance can be expressed as [28]where scattered voltage can be obtained using Generalized Telegraphers’ equations in the frequency domain with subsequent Inverse Fourier Transform (IFFT) [28].
Figure 14 shows the results obtained from simulations [27]. Obviously, high convergence rates are achieved for current mean since 5 × 5 points (i.e., 5 points for each RV) almost overlap data involving 7 × 7 SC points. Current variance provides important information regarding the data dispersion around mean values (via standard deviation in Figure 14). Therefore, the standard deviation is around 20 mA at = 100 μs.
A multivariate test case for the calculation of transient impedance is presented in Figure 15. The influence of variability of soil conductivity , lightning front time , and lightning time-to-half is investigated where these three input variables are modeled as random, each prescribed with uniform distribution as follows: (in mS/m), (in μs), and (in μs) [28]. The assessment of the mean value and variance using the SC full-tensor model with collocation points is given.
The conductivity of the medium has the highest impact during the entire simulation period, while front time has higher impact at the early time period and time-to-half is more dominant in the steady state.
4. Conclusion
The paper reviews some applications of stochastic-deterministic modeling in various areas of computational electromagnetics (CEM), previously reported by the authors.
Having completed the deterministic modeling, simple stochastic collocation (SC) formalism is implemented to efficiently account for uncertainties and to determine confidence intervals in the set of obtained numerical results.
The expansion of the statistical output in terms of mean and variance over a polynomial basis via stochastic collocation (SC) is proved to be a powerful method for uncertainty propagation and related sensitivity analysis demonstrating robustness and efficiency and providing satisfactory convergence rate and nonintrusiveness nature of the approach.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.