# Stochastic Collocation Applications in Computational Electromagnetics

Authors: Dragan Poljak; Silvestar Šesnić; Mario Cvetković; Anna Šušnjara; Hrvoje Dodig; Sébastien Lalléchère; Khalil El Khamlichi Drissi
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 Section2 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 Section3, 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 parameterX. Within this framework, a given mapping through function f:R→R can be assumed and a random variable (RV) X given by an a priori probabilistic distribution (f 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 functionf over nth order Lagrange polynomial basis can be written as follows:(1)fx=∑i=0nfiLix,where Lagrange polynomials are given by(2)Lix=∏k=0,k≠inx-xkxi-xk.And one has the following property:(3)Lixi=δij,where δij 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 pointsxi 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] yields(4)∫-∞+∞gxpxdx=∑j=0nωjgxj,where p is the probability density of random variable X, while g represents a function with a sufficient regularity and the coefficients ωj are collocation weights chosen in a way to ensure an exact quadrature rule for polynomials with a degree lower or equal to 2n+1.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) off(X) is by definition(5)EfX=∫-∞+∞fxpxdx,where X is a given RV.Next step is to replacef(x) in terms of Lagrange polynomials and expectation of f(X) can be written as follows:(6)EfX=∫-∞+∞∑i=0nfiLixpxdx.The Gauss quadrature rule yields(7)∫-∞+∞Lixpxdx=∑j=0nωjLixj.Taking into account the property of Lagrange polynomials (3) one obtains(8)∫-∞+∞Lixpxdx=ωi.Combining (6)–(8) leads to an expression for mean of f(X):(9)EfX=∑i=0nωifi.Furthermore, the same approach is used to evaluate variance off(X), which is by definition given by(10)varfX=EfX2-EfX2.It can be written as(11)varfX=∫-∞+∞fx2pxdx-EfX2.Again, utilizing the expansion of functionf over Lagrange polynomial basis (1) yields(12)varfX=∫-∞+∞∑i=0nfiLix∑j=0nfjLjxpxdx-EfX2which can be also written as(13)varfX=∑i=0n∑j=0nfifj∫-∞+∞LixLjxpxdx-EfX2.Furthermore, integral from the right-hand side of (13) according to Gauss quadrature rule (7) becomes(14)∫-∞+∞LixLjxpx=∑k=0nωkLixkLjxk=∑k=0nωkδikδjk.Finally, inserting (14) into (13), one obtains(15)varfX=∑i=0n∑j=0n∑k=0nωkδikδjk-EfX2 and due to Kronecker property nonzero terms are only with i=j=k leading to(16)varfX=∑k=0nωkfk2-EfX2.Now the statistical moments around a considered outputI can be written as Iu, where u stands for a given order of statistical moment; for example, the mean of f(X) can then be written as(17)EfX=I1.General expression is given as follows:(18)Iu≈∑i=0nωiΦuxi,where Φu represents a mapping pertaining to a certain statistical moment u, that is, with Φuxi such as(19)Φ1xi=fxi,Φ2xi=fxi-EfX2.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 vectorX=X1,X2,…,XNT and uth statistical moment of output I can be written as [28](20)IuX≈∑i1=0n1⋯∑iN=0nNωi1⋯iNΦuX,where ωsi1⋯iN stands for weight of the expansion for random vector X of size N. 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, Section3.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 lengthL and radius a and is located at height h above an interface as depicted in Figure 1.Figure 1
Center-fed GPR dipole antenna.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 currentI(x,t) along the dipole antenna excited by an equivalent voltage generator is obtained by solving the Hallen integral equation [31]:(21)∫0LIx′,t-Ra/c4πdx′-∫-∞t∫0Lrθ,τIx′,t-Ra∗/c-τ4πRa∗dx′dτ=12Z0Vgxg,t-x-xgc+F0t-xc+FLt-L-xc,where c is velocity of light, Z0 is the free-space wave impedance, and Ra and Ra∗ are the distances between observation point x and source point x′ on real and image wire, respectively. Time signals F0 and FL 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) r(θ,t) [31]:(22)rθ,τ=Kδt+4β1-β2e-αtt∑n=1∞-1n+1nKnInαt,where(23)K=1-β1+β;β=εr-sin2⁡θεrcos⁡⁡θ;α=σ2ε;θ=arctgx-x′2h;δ ( t ) is the Dirac impulse; and In(t) is the nth 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](24)Extrr,t=μ04π∫-∞t∫0LΓtrMITτ∂Ix′,t-R′′/v-τ∂τe-1/τgR′′/vR′′dx′dτ,where v is velocity of wave propagation in the lower medium and R′′ is the distance from the dipole antenna to the observation point located in the lower medium:(25)R′′=x-x′2+z+h2.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]:(26)Γtrt=τ3τ2δt+1τ22-τ3τ2e-t/τ2,where(27)τ2=εr+1σε0,τ3=2εσ.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 lengthL = 1 m and radius a = 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:(28)Vt=V0e-g2t-t02with parameters V0 = 1 V, g = 1.5∗109 s−1, t0 = 1.43 ns.Figure2 shows the mean of the electric field transmitted in the soil for the following random variables: RV1=h (antenna height above a soil), RV2=εr (soil permittivity), and RV3=d (field observation point). First, one random variable at a time is considered, as follows:(i)
RV1: central location, h0=0.3 m; uniformly distributed: U[0.05;0.55] (in meters);(ii)
RV2: initial value, εr0=15; uniformly distributed: U[2;28];(iii)
RV3: initial value, d0 = 1 m; uniformly distributed: U[0.5;1.5] (in meters).Figure 2
Mean of theE-field inside soil, L = 1 m, a = 6.74 mm. Figure reproduced from Poljak et al. (2017).Significant discrepancies appear between different scenarios with one random variable at a time, that is, betweenRV1, RV2, and RV3 models. The influence of each RV on the output E is computed using the following relation:(29)InfluenceEi=varE∣RVivarE∣RVi=1⋯n.Figure3 shows the influence of each of three different RVs on the output E.Figure 3
Influence of each RV on the outputE.Figure reproduced from Poljak et al. (2017).Ranking random parameters from the most influential to the least influential (fromE-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: εr0~U[7,30] and εr0~U[14,18] for wet and average sand, respectively. The other RVs are modeled as σ0~U[0.1,9.9] mS/m; h0~U[12,18] cm; d0~U[0.93,1.07] 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).Figure4 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 4
Electric field transmitted into the soil versus time for wet sand;Etr is mean value and std is standard deviation; only εr is random. Figure reproduced from Susnjara et al. (2017).Figure5 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.Figure 5
The impact factorIi(E) for each random input variable for the case of wet sand (a) and average type of sand (b). Figure reproduced from Susnjara et al. (2017).
(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 energyW absorbed by the unit body mass m:(30)SAR=dPdm=σ2ρE2, where P is the dissipated power, E 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]:(31)∇·k∇T+ρbcbwTa-T+Qm+ρ·SAR=0,where k stands for the tissue thermal conductivity, ρb is the blood mass density, cb is the specific heat capacity of blood, w is the blood perfusion rate, Ta is the arterial temperature, Qm is the heat source due to metabolic processes, and term ρ·SAR 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 is(32)q=HTs-Ta,where q is the heat flux density:(33)q=-λ∂T∂nwhile H, Ts, and Ta 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]:(34)jωμn∬SJ→r→′Gnr→,r→′dS′-jωεn∬S∇′SJ→r→′∇Gnr→,r→′dS′+∬SM→r→′×∇′Gnr→,r→′dS′=E→inc,n=10,n=2,where J→ and M→ are the unknown equivalent electric and magnetic current densities, respectively, kn is the wave number of a medium n, and Gn is the interior/exterior Green function for the homogeneous medium [32] given by(35)Gnr→,r→′=e-jknR4πR;R=r→-r→′,where R is the distance from the source to the observation point, respectively, while n denotes the index of an interior/exterior domain; that is, n=1,2.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]:(36)αE→=E→i+∮∂V′n→×∇′×E→GdS′+∮∂V′n→′×E→×∇′G+n→·E→GdS′,where Ei is the incident electric field.The interior domain is governed by the vector Helmholtz equation [25]:(37)∇×1kB∇×E→-kAE→=0,where subscripts A and B 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 atf = 900 MHz and f = 1800 MHz, respectively, for the case of vertical polarization. The power density of the incident plane wave is P=5 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 Table1.Table 1
The brain electrical parameters.
f = 900 MHz f = 1800 MHz ε r 46 84 σ [S/m] 0.8 1.2Figure6 shows the results for the SAR and the resulting temperature increase obtained in the brain model at 900 MHz [24].Figure 6
Distribution of SAR and the temperature increase due to 900 MHz vertically polarized incident plane wave. Figure reproduced from Cvetkovic et al. (2016).EachRVk(k=1,…,5), 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].Figure7 shows the SC convergence with 3, 5, and 7 points for the SAR values calculations and different random variables (RVk,k=1,…,5).Figure 7
Mean of maximum SAR. Figure reproduced from Cvetkovic et al. (2016).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.Figures8 and 9 are related to the results obtained from stochastic dosimetry modeling of the eye.Figure 8
Relative gap between SC approximations: (a)n=2 (SC3) and n=4 (SC5), (b) n=4 (SC5) and n=6 (SC7), regarding mean and variance from the real part of the induced electric field. Figure reproduced from Dodig et al. (2014).
(a) (b)Figure 9
Variance inside the eye model induced by plane wave with power density 10 W/m2 at f=6 GHz and SC approximation order n=2, results at sagittal cross section. Stochastic variations around vitreous body (a) conductivity, (b) relative permittivity, cornea (c) conductivity, and (d) relative permittivity. Figure reproduced from Dodig et al. (2014).
(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.Table 2
The electrical parameters of selected eye tissues.
Vitreous body Cornea ε r 64.5 ± 5 46.2 ± 5 σ [S/m] 7.01 ± 0.6 5.91 ± 0.6Figure8 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.Figure9 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 Figures9(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 lengthL and radius a is buried in a lossy medium at depth d, excited by the plane wave as shown in Figure 10.Figure 10
A horizontal thin wire scatterer buried in a lossy medium excited via plane wave.A transient response of a horizontal buried wire is governed by the following equation of the Pocklington type [26]:(38)με∂∂t+μσExtrt=-∂2∂x2-μσ∂∂t-με∂2∂t2·μ4π∫0Lix′,t-Rve-1/τgR/vRdx′-μ4π∫0t∫0LΓrefMITτix′,t-R∗v-τe-1/τgR∗/vR∗dx′dτ,where ix′,t-R/v is the unknown space-time dependent current, Extr is the tangential component of the transmitted field, and ΓrefMIT 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 is(39)R=x-x′2+a2while the distance from the source point on the image wire to the observation point on the original wire, according to the image theory, is(40)R∗=x-x′2+4d2.Time constant and propagation velocity in the lossy medium are given with(41)τg=2εσ,v=1με.The influence of the earth-air interface is taken into account via the reflection coefficient arising from the MIT and is given with [34](42)ΓrefMITt=-τ1τ2δt+1τ21-τ1τ2e-t/τ2,where the corresponding time constants are(43)τ1=ε0εr-1σ,τ2=ε0εr+1σ.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σ: σ0=X10=5 mS/m; uniformly distributed: U[1;9] (in mS/m);(ii)
lengthL of wire: L0=X20=10 m; uniformly distributed: U[9.5;10.5] (in meters);(iii)
burial depthd:  d0=X30=4 m; uniformly distributed: U[2.5;5.5] (in meters).Without loss of generality, the problem can be addressed following different assumptions about the statistical distribution laws.Figure11 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 11
Currents at the center of the wire (“full-tensor” model). Figure reproduced from Sesnic et al. (2015).Figure12 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 d (RV3). Nearly zero levels of the current (mean and variance) below 0.03 μs involve instability of the convergence criterion (and positive SC gaps).Figure 12
Relative gap with increasing SC orders (RV3). Figure reproduced from Sesnic et al. (2015).


### 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 lengthL and radius a, buried in a lossy medium of conductivity σ and relative permittivity ε, at depth d, is shown in Figure 13.Figure 13
A horizontal thin wire scatterer buried in a lossy medium excited via equivalent current source.Thus, the current distribution along the wire is governed by the space-time homogeneous Pocklington integrodifferential equation [35]:(44)∂2∂x2-μσ∂∂t-με∂2∂t2·μ4π∫0Lix′,t-Rve-1/τgR/vRdx′-μ4π∫0t∫0LΓrefMITτix′,t-R∗v-τe-1/τgR∗/vR∗dx′dτ.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](45)zt=v0,ti0,t,where scattered voltage can be obtained using Generalized Telegraphers’ equations in the frequency domain with subsequent Inverse Fourier Transform (IFFT) [28].Figure14 shows the results obtained from 32+52+72=83 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 t = 100 μs.Figure 14
Mean and standard deviation of the current at the center of the electrode in function of time. Figure reproduced from Sesnic et al. (2016).A multivariate test case for the calculation of transient impedance is presented in Figure15. The influence of variability of soil conductivity σ, lightning front time TF, and lightning time-to-half TH is investigated where these three input variables are modeled as random, each prescribed with uniform distribution as follows: σ~U(1,10) (in mS/m), TF~U(0.4,4) (in μs), and TH~U(50,70) (in μs) [28]. The assessment of the mean value and variance using the SC full-tensor model with 7×7×7 collocation points is given.Figure 15
Mean value ± standard deviation of transient impedance for multivariate case. Figure reproduced from Sesnic et al. (2017).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.


## 3.1. Transient Field in the Soil Generated by GPR Dipole Antenna

The dipole is characterized by a lengthL and radius a and is located at height h above an interface as depicted in Figure 1.Figure 1
Center-fed GPR dipole antenna.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 currentI(x,t) along the dipole antenna excited by an equivalent voltage generator is obtained by solving the Hallen integral equation [31]:(21)∫0LIx′,t-Ra/c4πdx′-∫-∞t∫0Lrθ,τIx′,t-Ra∗/c-τ4πRa∗dx′dτ=12Z0Vgxg,t-x-xgc+F0t-xc+FLt-L-xc,where c is velocity of light, Z0 is the free-space wave impedance, and Ra and Ra∗ are the distances between observation point x and source point x′ on real and image wire, respectively. Time signals F0 and FL 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) r(θ,t) [31]:(22)rθ,τ=Kδt+4β1-β2e-αtt∑n=1∞-1n+1nKnInαt,where(23)K=1-β1+β;β=εr-sin2⁡θεrcos⁡⁡θ;α=σ2ε;θ=arctgx-x′2h;δ ( t ) is the Dirac impulse; and In(t) is the nth 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](24)Extrr,t=μ04π∫-∞t∫0LΓtrMITτ∂Ix′,t-R′′/v-τ∂τe-1/τgR′′/vR′′dx′dτ,where v is velocity of wave propagation in the lower medium and R′′ is the distance from the dipole antenna to the observation point located in the lower medium:(25)R′′=x-x′2+z+h2.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]:(26)Γtrt=τ3τ2δt+1τ22-τ3τ2e-t/τ2,where(27)τ2=εr+1σε0,τ3=2εσ.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 lengthL = 1 m and radius a = 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:(28)Vt=V0e-g2t-t02with parameters V0 = 1 V, g = 1.5∗109 s−1, t0 = 1.43 ns.Figure2 shows the mean of the electric field transmitted in the soil for the following random variables: RV1=h (antenna height above a soil), RV2=εr (soil permittivity), and RV3=d (field observation point). First, one random variable at a time is considered, as follows:(i)
RV1: central location, h0=0.3 m; uniformly distributed: U[0.05;0.55] (in meters);(ii)
RV2: initial value, εr0=15; uniformly distributed: U[2;28];(iii)
RV3: initial value, d0 = 1 m; uniformly distributed: U[0.5;1.5] (in meters).Figure 2
Mean of theE-field inside soil, L = 1 m, a = 6.74 mm. Figure reproduced from Poljak et al. (2017).Significant discrepancies appear between different scenarios with one random variable at a time, that is, betweenRV1, RV2, and RV3 models. The influence of each RV on the output E is computed using the following relation:(29)InfluenceEi=varE∣RVivarE∣RVi=1⋯n.Figure3 shows the influence of each of three different RVs on the output E.Figure 3
Influence of each RV on the outputE.Figure reproduced from Poljak et al. (2017).Ranking random parameters from the most influential to the least influential (fromE-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: εr0~U[7,30] and εr0~U[14,18] for wet and average sand, respectively. The other RVs are modeled as σ0~U[0.1,9.9] mS/m; h0~U[12,18] cm; d0~U[0.93,1.07] 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).Figure4 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 4
Electric field transmitted into the soil versus time for wet sand;Etr is mean value and std is standard deviation; only εr is random. Figure reproduced from Susnjara et al. (2017).Figure5 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.Figure 5
The impact factorIi(E) for each random input variable for the case of wet sand (a) and average type of sand (b). Figure reproduced from Susnjara et al. (2017).
(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 energyW absorbed by the unit body mass m:(30)SAR=dPdm=σ2ρE2, where P is the dissipated power, E 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]:(31)∇·k∇T+ρbcbwTa-T+Qm+ρ·SAR=0,where k stands for the tissue thermal conductivity, ρb is the blood mass density, cb is the specific heat capacity of blood, w is the blood perfusion rate, Ta is the arterial temperature, Qm is the heat source due to metabolic processes, and term ρ·SAR 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 is(32)q=HTs-Ta,where q is the heat flux density:(33)q=-λ∂T∂nwhile H, Ts, and Ta 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]:(34)jωμn∬SJ→r→′Gnr→,r→′dS′-jωεn∬S∇′SJ→r→′∇Gnr→,r→′dS′+∬SM→r→′×∇′Gnr→,r→′dS′=E→inc,n=10,n=2,where J→ and M→ are the unknown equivalent electric and magnetic current densities, respectively, kn is the wave number of a medium n, and Gn is the interior/exterior Green function for the homogeneous medium [32] given by(35)Gnr→,r→′=e-jknR4πR;R=r→-r→′,where R is the distance from the source to the observation point, respectively, while n denotes the index of an interior/exterior domain; that is, n=1,2.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]:(36)αE→=E→i+∮∂V′n→×∇′×E→GdS′+∮∂V′n→′×E→×∇′G+n→·E→GdS′,where Ei is the incident electric field.The interior domain is governed by the vector Helmholtz equation [25]:(37)∇×1kB∇×E→-kAE→=0,where subscripts A and B 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 atf = 900 MHz and f = 1800 MHz, respectively, for the case of vertical polarization. The power density of the incident plane wave is P=5 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 Table1.Table 1
The brain electrical parameters.
f = 900 MHz f = 1800 MHz ε r 46 84 σ [S/m] 0.8 1.2Figure6 shows the results for the SAR and the resulting temperature increase obtained in the brain model at 900 MHz [24].Figure 6
Distribution of SAR and the temperature increase due to 900 MHz vertically polarized incident plane wave. Figure reproduced from Cvetkovic et al. (2016).EachRVk(k=1,…,5), 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].Figure7 shows the SC convergence with 3, 5, and 7 points for the SAR values calculations and different random variables (RVk,k=1,…,5).Figure 7
Mean of maximum SAR. Figure reproduced from Cvetkovic et al. (2016).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.Figures8 and 9 are related to the results obtained from stochastic dosimetry modeling of the eye.Figure 8
Relative gap between SC approximations: (a)n=2 (SC3) and n=4 (SC5), (b) n=4 (SC5) and n=6 (SC7), regarding mean and variance from the real part of the induced electric field. Figure reproduced from Dodig et al. (2014).
(a) (b)Figure 9
Variance inside the eye model induced by plane wave with power density 10 W/m2 at f=6 GHz and SC approximation order n=2, results at sagittal cross section. Stochastic variations around vitreous body (a) conductivity, (b) relative permittivity, cornea (c) conductivity, and (d) relative permittivity. Figure reproduced from Dodig et al. (2014).
(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.Table 2
The electrical parameters of selected eye tissues.
Vitreous body Cornea ε r 64.5 ± 5 46.2 ± 5 σ [S/m] 7.01 ± 0.6 5.91 ± 0.6Figure8 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.Figure9 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 Figures9(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 lengthL and radius a is buried in a lossy medium at depth d, excited by the plane wave as shown in Figure 10.Figure 10
A horizontal thin wire scatterer buried in a lossy medium excited via plane wave.A transient response of a horizontal buried wire is governed by the following equation of the Pocklington type [26]:(38)με∂∂t+μσExtrt=-∂2∂x2-μσ∂∂t-με∂2∂t2·μ4π∫0Lix′,t-Rve-1/τgR/vRdx′-μ4π∫0t∫0LΓrefMITτix′,t-R∗v-τe-1/τgR∗/vR∗dx′dτ,where ix′,t-R/v is the unknown space-time dependent current, Extr is the tangential component of the transmitted field, and ΓrefMIT 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 is(39)R=x-x′2+a2while the distance from the source point on the image wire to the observation point on the original wire, according to the image theory, is(40)R∗=x-x′2+4d2.Time constant and propagation velocity in the lossy medium are given with(41)τg=2εσ,v=1με.The influence of the earth-air interface is taken into account via the reflection coefficient arising from the MIT and is given with [34](42)ΓrefMITt=-τ1τ2δt+1τ21-τ1τ2e-t/τ2,where the corresponding time constants are(43)τ1=ε0εr-1σ,τ2=ε0εr+1σ.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σ: σ0=X10=5 mS/m; uniformly distributed: U[1;9] (in mS/m);(ii)
lengthL of wire: L0=X20=10 m; uniformly distributed: U[9.5;10.5] (in meters);(iii)
burial depthd:  d0=X30=4 m; uniformly distributed: U[2.5;5.5] (in meters).Without loss of generality, the problem can be addressed following different assumptions about the statistical distribution laws.Figure11 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 11
Currents at the center of the wire (“full-tensor” model). Figure reproduced from Sesnic et al. (2015).Figure12 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 d (RV3). Nearly zero levels of the current (mean and variance) below 0.03 μs involve instability of the convergence criterion (and positive SC gaps).Figure 12
Relative gap with increasing SC orders (RV3). Figure reproduced from Sesnic et al. (2015).


## 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 lengthL and radius a, buried in a lossy medium of conductivity σ and relative permittivity ε, at depth d, is shown in Figure 13.Figure 13
A horizontal thin wire scatterer buried in a lossy medium excited via equivalent current source.Thus, the current distribution along the wire is governed by the space-time homogeneous Pocklington integrodifferential equation [35]:(44)∂2∂x2-μσ∂∂t-με∂2∂t2·μ4π∫0Lix′,t-Rve-1/τgR/vRdx′-μ4π∫0t∫0LΓrefMITτix′,t-R∗v-τe-1/τgR∗/vR∗dx′dτ.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](45)zt=v0,ti0,t,where scattered voltage can be obtained using Generalized Telegraphers’ equations in the frequency domain with subsequent Inverse Fourier Transform (IFFT) [28].Figure14 shows the results obtained from 32+52+72=83 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 t = 100 μs.Figure 14
Mean and standard deviation of the current at the center of the electrode in function of time. Figure reproduced from Sesnic et al. (2016).A multivariate test case for the calculation of transient impedance is presented in Figure15. The influence of variability of soil conductivity σ, lightning front time TF, and lightning time-to-half TH is investigated where these three input variables are modeled as random, each prescribed with uniform distribution as follows: σ~U(1,10) (in mS/m), TF~U(0.4,4) (in μs), and TH~U(50,70) (in μs) [28]. The assessment of the mean value and variance using the SC full-tensor model with 7×7×7 collocation points is given.Figure 15
Mean value ± standard deviation of transient impedance for multivariate case. Figure reproduced from Sesnic et al. (2017).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.



---

*Source: 1917439-2018-05-14.xml*