Maryvonne Giron

Céline Labrude

Frédéric Desprez

Keywords: equivalent liquid, 3, 4 SAR calculation at 890 MHz using different normalizations

Three years ago, the China Scholarship Council (CSC) awarded me a three-years scholarship to work abroad as Ph.D student. At that time, my Ph.D topic was entirely new and quite mysterious to me, since my Master topic was completely different. And as a foreigner who started to learn French at the age of 25 and only for three months, communicating with people in French was difficult. Nevertheless, Professor Vikass Monebhurrun kindly received me in his research group, where he patiently taught me electromagnetics and in particular the human exposure to electromagnetic fields. I would like to thank him very sincerely at first, for his scientific advice at that time and since then. Then I would like to thank Dr. Dominique Lesselier for his kind support during the three years of my Ph.D. Even though I am not his Ph.D student, and I have worked in the same office with his students, Dr. Dominique Lesselier has treated me as his student and helped me sincerely. For research, he gave me some advice and provided me papers which related to my topic. Two years ago, since I had difficulty in communicating with people in French, he helped me to apply to the temporary residence on France as a Ph.D student. Since there were a lots of documents to deal with, he spent much time to collect these documents for me. Here I would like to thank him very much.

I much appreciate the help of the people who work in the laboratory,

Introduction et contexte de recherche

La procédure de mesure de la conformité des téléphones mobiles en termes de débit d'absorption spécifique (DAS, ou Specific Absorption Rate, SAR) [1-6] est décrite par les normes internationales. Puisque le test de conformité DAS d'un téléphone mobile est coûteux et prend beaucoup de temps, la disponibilité de logiciels de simulation électromagnétique fonctionnant sur des ordinateurs standard offre des perspectives intéressantes pour le calcul de DAS des téléphones mobiles.

Une procédure standardisée pour le calcul de DAS des téléphones mobiles par l'utilsation de la méthode des différences finies dans le domaine temporel (FDTD) [7] est actuellement développée au sein d'un groupe de travail de l'ICES/IEEE/TC34 tel que considéré en [8,9]. Un benchmark de Conception Assistée par Ordinateur (CAO) de modèle de téléphone mobile a été récemment proposé pour une comparaison internationale inter-laboratoires d'outils de simulation électromagnétiques [10]. Bien qu'il ait été obtenu un bon accord global entre les résultats de neuf laboratoires, un écart maximal d'environ 20 % a été observé pour ce qui concerne le maximum du DAS moyenne sur 10g.

Puisque le même modèle de téléphone mobile et les mêmes paramètres de matériaux ont été utilisés par les différents laboratoires participants, parfois même les mêmes logiciels utilisés, les différences dans les densités de maillage choisies et le positionnement par rapport au fantôme de tête pourraient expliquer les écarts observés. En résumé, afin d'obtenir des résultats faisant sens, il est important d'évaluer l'incertitude des calculs de DAS. Mais la quantification des incertitudes (Uncertainty Quantification ou UQ) dans le calcul de DAS en utilisant un modèle de téléphone mobile réaliste est une tâche difficile.

La Commission internationale sur la protection contre les rayonnements non-ionisants (International Commission on Non-Ionizing Radiation Protection ou ICNIRP) [4] a analysé en détail l'absorption de l'énergie au sein du corps humain causée par des champs électromagnétiques. Il est en particulier expliqué que l'exposition aux champs basses fréquences électriques et magnétiques se traduit normalement par une absorption d'énergie négligeable et point de hausse de température mesurable dans le corps. Toutefois, l'exposition aux champs électromagnétiques à des fréquences supérieures à 100 kHz peut, elle, conduire à une absorption importante de l'énergie et la température augmenter en conséquence.

En général, l'exposition à un champ électromagnétique uniforme (type onde plane) conduit au dépôt et à une distribution de l'énergie hautement non-uniformes dans le ix corps, phénomène qui doit être évalué tant par une mesure dosimétrique que par le calcul. En ce qui concerne l'absorption d'énergie par l'organisme humain, les champs électromagnétiques peuvent être considérés via quatre plages de fréquences [4] :

• fréquences de l'ordre de 100 kHz à moins d'environ 20 MHz, pour lesquelles l'absorption dans le tronc décroît rapidement avec la diminution de la fréquence, tandis qu'une absorption significative peut se produire dans le cou et les jambes ;

• fréquences d'environ 20 MHz à 300 MHz, plage au sein de laquelle une absorption relativement élevée peut se produire dans tout le corps, pouvant peut-être atteindre des valeurs encore plus élevées si une partie du corps (par exemple, la tête) est telle que des résonances apparaissent ;

• fréquences d'environ 300 MHz à plusieurs GHz, plage au sein de laquelle une absorption significative, non localement uniforme, se produit ;

• fréquences supérieures à 10 GHz, qui conduisent à ce que l'absorption d'énergie s'effectue principalement à la surface du corps.

Le DAS est proportionnel au carré de l'intensité du champ électrique dans le tissu interne. Une répartition moyenne et une distribution de DAS peuvent être calculées ou estimées à partir de mesures de laboratoire. Les valeurs de DAS dépendent des facteurs suivants :

• paramètres du champ incident, à savoir la fréquence, l'intensité, la polarisation et la configuration objet-source (notions de champ proche ou lointain) ;

• caractéristiques du corps exposé, à savoir sa taille et sa géométrie tant interne qu'externe, et les propriétés diélectriques des différents tissus ;

• effets de masse et effets de réflecteur d'autres objets présents dans le champ proche du corps exposé.

L'utilité de la modélisation numérique ainsi que des mesures de courant induit dans le corps et de l'intensité du champ au sein des tissus pour l'évaluation de l'exposition en champ proche a été prouvée pour téléphones mobiles, tours d'antennes de radiodiffusion, sources de communication à bord des navires, systèmes de chauffage de diélectriques, etc. L'ICNIRP a ensuite établi des directives afin de limiter l'exposition aux champs électromagnétiques et a défini des contraintes de base qui spécifient un DAS ne devant pas être dépassé et les niveaux de référence en termes de champ électrique, champ x Résumé magnétique, ou densité de puissance. Les niveaux de référence ont été dérivés des contraintes de base, il y a plusieurs années, à partir d'un modèle humain très simplifié. Depuis lors, un certain nombre de modèles humains réalistes a été développé et utilisé en dosimétrie.

En [11], un bref historique de la modélisation et les détails de la quantification d'incertitude en modélisation sont proposés. Ci-après des éléments pertinents en sont donnés.

Le modèle d'un phénomène physique est défini comme étant une représentation abstraite idéalisée de ce phénomène par un ensemble d'équations dont la solution reproduit les observations expérimentales. La fonction d'un modèle est que les ingénieurs obtiennent les lois de la physique, par exemple, gravitation, de l'optique, de l'électromagnétisme, etc., à partir des prévisions du modèle. Si modélisation, il y a erreur, erreur numérique et incertitude existant dans la modélisation. L'erreur de modèle vient d'un manque de connaissance de la physique vraie sous-jacente. Elle dépend de à quel niveau de précision un modèle mathématique décrit le système vrai pour une situation de la vie réelle, compte tenu du fait que les modèles sont presque toujours des approximations de la réalité seulement.

Les modèles mathématiques sont en règle générale trop compliqués pour être traités exactement. Certains algorithmes peuvent être utilisés pour se rapprocher de la solution, ce qui, cependant, introduit des erreurs numériques ; par exemple, la méthode des éléments finis (Finite Element Method ou FEM) [12] ou la méthode des différences finies (Finite Difference Method ou FDM) [13] peuvent être utilisées pour se rapprocher de la solution d'une équation aux dérivées partielles. La définition de l'incertitude est une carence potentielle quelle que soit la phase ou l'activité du processus de modélisation, par un manque de connaissances. L'incertitude peut elle-même être classée en deux catégories [14] : L'incertitude de "Type A" ou incertitude aléatoire résulte du comportement stochastique ou du comportement aléatoire dans le système. Elle n'existe pas dans des simulations de calcul déterministes.

L'incertitude de "Type B" ou incertitude épistémique est l'incertitude qui se pose en raison du manque de connaissances. Ce type d'incertitude peut être réduit par l'acquisition de plus de connaissances. Elle peut être en outre classée en deux groupes.

Le premier groupe d'incertitude épistémique est l'incertitude sur la façon dont le modèle mathématique représente le vrai comportement du système physique. Cette incertitude de forme du modèle est très difficile à déterminer, mais peut être réduite en vérifiant le modèle par le biais de mesures physiques.

xi

Le second groupe est l'incertitude qui apparaît en raison d'un manque de données de paramètres d'entrée précis. Pour un système, les incertitudes dans les paramètres d'entrée vont induire les incertitudes dans les sorties. Ce type d'incertitude est souvent qualifiée en incertitude paramétrique.

Méthodologies et expérimentations numériques

La méthode de Monte Carlo (Monte Carlo Method ou MCM) est la méthode la plus reconnue pour estimer l'incertitude dans diverses simulations de compatibilité électromagnétique (CEM). L'inconvénient majeur de la MCM est qu'elle ne permet que des convergences lentes. Il est ainsi impossible de la mettre en oeuvre dans les calculs de DAS en raison de cet inconvénient. En effet, un calcul typique de DAS d'un téléphone mobile placé contre une tête fantôme de Mannequin Antropomorphe Spécifique (Specific Anthropomorphic Mannequin ou SAM) nécessite environ une demi-heure ou même plusieurs heures en fonction de la machine de calcul utilisée pour la simulation.

Cela est évidemment très peu pratique si l'on doit effectuer plusieurs milliers de calculs de DAS, et certaines méthodes alternatives doivent donc être proposées pour la quantification d'incertitude (UQ) dans le calcul de DAS. Dans le chapitre correspondant (le chapitre 2) de la thèse, l'une des méthodes efficaces de calcul de UQ est proposée à cette fin.

Cette méthode est appelée Transformation sans parfum (Unscented Transformation ou UT). En substance, UT est assez similaire à MCM puisque les deux méthodes sont utilisées afin d'obtenir une approximation de la fonction de distribution de probabilité plutôt que la fonction non-linéaire elle-même, qui est donc considérée comme une boîte noire. UT est une méthode non intrusive d'évaluation qui est proposée dans [15] et ensuite utilisée pour le calcul d'UQ dans le domaine de la CEM [START_REF] Menegaz | A new smallest sigma set for the Unscented Transform and its applications on SLAM[END_REF]17,18]. L'idée de l'UT est présentée en détail dans ledit chapitre, et les résultats de l'application de UT et MCM pour l'UQ dans les calculs sont comparés et analysés.

Deux autres méthodes UQ non intrusives sont introduites ensuite dans le chapitre 3. Ces deux méthodes sont la méthode de collocation stochastique (Stochastic Collocation ou SC) et celle des polynômes de chaos non-intrusifs (Non-Intrusive Polynomial Chaos ou NIPC).

Les méthodes SC et PC sont des méthodes en fait très similaires parce que toutes les deux ont la capacité de produire des représentations fonctionnelles de la variabilité stochastique. La méthode SC construit des fonctions d'interpolation pour les coefficients connus, et exige l'utilisation d'ensembles de points de collocation structurés issus de xii Résumé produits tensoriels ou grilles parcimonieuses. La méthode PC, quant à elle, estime des coefficients pour des fonctions de base connues de polynômes orthogonaux en se fondant sur une série d'évaluations de la fonction de réponse, en utilisant échantillonnage, régression linéaire, intégration de produits tensoriels, ou des approches de Smolyak sur grilles parcimonieuses.

Les performances des méthodes UT, SC et NIPC pour la quantification d'incertitude dans le calcul de DAS en utilisant un modèle de téléphone mobile commercial disponible par ailleurs sont comparées et analysées. Pour quantifier l'importance relative des entrées incertaines dans le calcul de DAS, l'une des méthodes est l'analyse de sensibilité (Sensitivity Analysis ou SA) basée sur la variance, en lien à l'expansion en polynômes de chaos, qui est utilisée dans ce but.

Toutes les méthodes d'UQ y compris les méthodes UT, SC et NIPC possèdent un inconvénient commun en raison du fait qu'il est difficile de les utiliser avec de multiples entrées incertaines dans le calcul. Les formulations de méthodes UQ multidimensionnelles sont complexes, et le nombre des simulations requises par l'UQ augmente avec le nombre de paramètres d'entrée. Et dans le calcul de DAS, il y a généralement plusieurs variables aléatoires en raison de l'absence de connaissance précise de la matière ou de la géométrie des composants du téléphone mobile étant modélisées. Pour résoudre ce problème, des incertitudes combinées sont étudiées pour être utilisées dans le calcul de DAS.

Avant que la méthode NIPC ne soit proposée et étudiée, une approche en polynômes de chaos généralisé (generalized Polynomial Chaos ou gPC) est utilisée de manière intrusive, pour UQ au premier abord. Puisque le développement gPC est appliqué dans le code de simulation tel que le code FDTD, les méthodes sont dites intrusives.

A ce sujet, en [19] est discuté un travail impliquant la mise en oeuvre de l'expansion gPC dans une solution de Galerkin discontinu d'ordre supérieur des équations de Maxwell, observant que la méthode gPC peut quantifier avec précision l'incertitude dans la sortie, tout en étant plus efficace que la MCM ; cependant seule l'incertitude de sortie pour un paramètre d'entrée incertain a été considérée. Puis, en [20], un nouvel algorithme d'application de gPC a été proposé et appliqué en FDTD afin de quantifier l'incertitude de la sortie d'une simulation CEM. Les équations de mise à jour du nouvel algorithme suivent une forme similaire à celle de celui de la FDTD conventionnelle. Les résultats de l'incertitude de sortie qui est induite par trois paramètres d'entrée incertains indépendants évalués par la nouvelle méthode pour un exemple simple de simulation CEM sont ainsi donnés. Une application est faite en [21] du nouvel algorithme pour mener à bien différents types de simulation CEM, dont xiii visant champ rétrodiffusé et section efficace radar en diffraction électromagnétique, pour analyser l'impact des incertitudes sur les résultats de simulation, et la méthode a prouvé son efficacité.

Dans le chaptre 4, le nouvel algorithme de l'application de la méthode gPC dans la simulation FDTD est introduit en détail. Avant l'introduction de la méthode, la connaissance nécessaire, telle que la formulation du système stochastique et la méthode FDTD, est fournie. Un exemple d'application de la méthode gPC dans une simulation 1D-FDTD est donné. Les résultats de la méthode sont comparées à ceux des méthodes MCM et NIPC.

Conclusion

Comme il est impossible de mettre en oeuvre la MCM dans le calcul de DAS en raison du fait qu'elle est trop coûteuse en calculs, investiguer des méthodes de calcul d'UQ efficaces devient indispensable. La principale contribution de cette thèse porte sur l'étude de certaines méthodes UQ non intrusives de calcul efficaces dans le calcul de DAS. L'objectif premier est un temps de calcul et/ou un nombre de simulations acceptables dans une évaluation d'incertitude dans le calcul de DAS, qui est similaire à celle du résultat obtenu par MCM.

Les méthodes non intrusives UQ appliquées dans cette thèse sont donc UT, SC et NIPC, et leur mise en oeuvre est aussi simple que celle de MCM, ne nécessitant que l'exécution répétitive d'un code déterministe à un ensemble de points prédéfinis dans l'espace stochastique. Ces méthodes sont satisfaisantes dans des conditions spécifiques. Elles donnent de bons résultats pour l'analyse d'incertitude, cependant, le nombre de simulations requises par UT augmente rapidement avec la dimension des vecteurs aléatoires.

Il convient donc de ne pas mettre en oeuvre l'UT quand il existe plusieurs paramètres d'entrée incertains dans le calcul. Lorsque les points de collecte de SC et les coefficients des polynômes de NIPC sont calculés par la règle de quadrature, il y a des erreurs existantes dans les résultats des deux méthodes en raison de cette règle de quadrature.

En utilisant la méthode de régression afin d'évaluer les coefficients polynomiaux dans NIPC, la mise en oeuvre d'un sur-échantillonnage est un bon choix pour améliorer les performances du procédé, même si l'efficacité de la méthode est diminuée en conséquence. L'erreur due à la troncature finie d'un développement infini existe en SC et NIPC.

Summary

This thesis focuses on parameter Uncertainty Quantification (UQ) in Specific Absorption Rate (SAR) calculation of mobile phone. The impact of uncertainty, e.g., lack of detailed knowledge about material electrical properties, system geometrical features, etc., in SAR calculation is quantified by three computationally efficient non-intrusive UQ methods: Unscented Transformation (UT), Stochastic Collocation (SC) and Non-Intrusive Polynomial Chaos (NIPC). They are called non-intrusive methods since the simulation process is simply considered as a black-box without changing the code of the simulation solver. Their performances for the cases of one and two random variables are analysed. In contrast with the traditional uncertainty analysis method (the Monte Carlo Method (MCM)), the computational cost becomes acceptable. To simplify the UQ procedure for the case of multiple uncertain inputs, it is demonstrated that uncertainties can be combined so as to evaluate the parameter uncertainty of the output. Combining uncertainties is an approach generally used in the field of measurement, in this thesis, it is used in SAR calculation for the complex situation. One of the necessary steps in the framework of uncertainty analysis is Sensitivity Analysis (SA) which aims at quantifying the relative importance of each uncertain input parameter with respect to the uncertainty of the output. Polynomial Chaos (PC) based Sobol' indices method the sensitivity indices of which are evaluated by PC expansion instead of by the MCM is used in SAR calculation. The results of the investigations are presented and discussed.

In order to ease the reading of the manuscript, elementary notions of SAR, modelling, uncertainty in modelling, probability theory, and SAR calculation in one of the simulation solvers are given in introduction chapter 1. Then the main content of this thesis is presented in chapter 2 and chapter 3. In chapter 4, another approach to use PC expansion is given, and it is used in the Finite-Difference Time-Domain (FDTD) code. Since the FDTD code in the simulation solver should be changed, it is the so-called intrusive PC expansion. Intrusive PC method has been investigated in detail by many researchers. In chapter 5, conclusions and an outline of future work are provided.

Table of contents

Résumé vii

References xv

List of figures xxiii

List of tables xxv

Uncertainty in SAR calculation

The measurement procedure of the Specific Absorption Rate (SAR) [1-7] compliance of mobile phones is described in international standards. Since the SAR compliance test of a mobile phone is costly and time-consuming, the availability of electromagnetic simulation software running on standard computers offers interesting perspectives.

A standardized procedure for the SAR calculation of mobile phones using the Finite-Difference Time-Domain (FDTD) method [8] is currently being developed by a working group of ICES/IEEE/TC34 [9,10]. A benchmark Computer-Aided Design (CAD) mobile phone model was recently proposed for an international inter-laboratory comparison of electromagnetic simulation tools [11]. Although there was a good overall agreement among the results of the nine laboratories, a maximum deviation of about 20% was observed in the peak 10g average SAR value. Yet, the same mobile phone model and material parameters have been used by the different laboratories, sometimes using the same software, so the differences in the applied mesh densities and positioning against the head phantom might explain the deviations observed.

To achieve meaningful results, it is important to evaluate the uncertainty of the SAR calculations. However, Uncertainty Quantification (UQ) in SAR calculation using a realistic mobile phone model is a rather challenging task. So, before entering into this topic in depth, the related knowledge is summarized in this chapter.

Specific absorption rate (SAR)

SAR is a measure of the rate at which energy is absorbed by the human body when exposed to a Radio Frequency (RF) electromagnetic field. SAR is usually averaged either over the whole body, or over a small sample volume (typically 1 Gram (g) or 10 Grams (g) of tissue). The value cited then is the maximum level measured in the body part studied over the stated volume or mass. SAR is defined as the power absorbed per mass of tissue and has units of watts per kilogram (W /kg) [4]

SAR = σE 2 ρ (1.1)
where σ is the conductivity of human tissue, ρ is its mass density and E is the modulus of the electric field induced in the tissue.

The International Commission on Non-Ionizing Radiation Protection (ICNIRP) [4] documents in detail the absorption of energy of the human body from electromagnetic fields. Experts in the Commission explain that exposure to low-frequency electric and magnetic fields normally results in negligible energy absorption and no measurable temperature rise in the body. However, exposure to electromagnetic fields at frequencies above about 100 kHz can lead to significant absorption of energy and temperature increases. In general, exposure to a uniform electromagnetic field results into a highly non-uniform deposition and distribution of energy within the body, to assess by dosimetric measurement and calculation.

As regards absorption of energy by the body, frequencies can be divided as follows:

• from about 100 kHz to less than about 20 MHz, for which absorption in the trunk decreases rapidly with decreasing frequency, and significant absorption may occur in the neck and legs;

• from about 20 MHz to 300 MHz, for which relatively high absorption can occur in the whole body, and even higher values if partial body (e.g., head) resonances are considered;

• from about 300 MHz to several GHz, for which significant local, non-uniform absorption occurs;

• frequencies above about 10 GHz, for which energy absorption primarily occurs at the body surface.

SAR is proportional to the square of the internal electric field strength in tissue. Average SAR and SAR distribution can be computed or estimated from laboratory measurements. Values of SAR mostly depend upon the incident field parameters (frequency, intensity, polarization, near or far field), the characteristics of the exposed body (size and internal and external geometry) and the dielectric properties of the various tissues, 1.1 Uncertainty in SAR calculation 3 and ground effects and reflector effects of other objects in the field near the exposed body.

The usefulness of numerical modelling calculations, as well as measurements of induced body current and tissue field strength, for assessment of near field exposures has been demonstrated for mobile telephones, walkie-talkies, broadcast towers, shipboard communication sources, dielectric heaters, in particular.

ICNIRP established guidelines for limiting electromagnetic field exposure and defined basic restrictions which specify a SAR not to be exceeded and reference levels in terms of electric field, magnetic field, or power density. The latter have been derived in effect from basic restrictions several years ago involving a very simplified human model. Since, many realistic human models have been developed for dosimetry studies [5].

Finite Integration Technique (FIT) and SAR calculation in CST

There are many commercial electromagnetic solvers, such as EMPIRE based on the FDTD method from IMST, EMPro based on the FDTD method from Agilent, and Microwave Studio based on the Finite Integral Technique (FIT) [START_REF] Clemens | Numerical simulation of coupled transient thermal and electromagnetic fields with the finite integration method[END_REF] from CST [11]. In this thesis, all SAR calculations are implemented by Microwave Studio from CST.

To remind, it is a general-purpose electromagnetic simulator based on FIT. FIT provides an universal spatial discretization scheme applicable to problems ranging from static to high frequency applications in time or frequency. The discretization in FIT is similar with the one in FDTD for homogeneous media. However, FIT transforms Maxwell's equations, in their integral form, into a linear system of equations. It treats interfaces in accurate manner, and is rather flexible in terms of geometry. In short, FIT discretizes the following integral form of these equations rather than the differential one:

E • ds = - ∂ B ∂ t • dA (1.2) D • dA = ρdV (1.3) H • ds = ∂ D ∂ t + J • dA (1.4) B • dA = 0 (1.5)
where E is the electric field, B is the magnetic field, dA denotes the differential vector element of surface area A, dV is the differential volume element, J is the electric current density, ds is the differential vector element of path length tangential to the path, D is the displacement field, H is the magnetizing field and ρ is the total charge density.

To solve these equations numerically, a finite calculation domain must be defined to enclose the considered application problem. CST Microwave Studio offers whole-body-averaged and local SAR values:

• Whole-body-averaged SAR: the value is obtained by dividing the total power absorbed in the human body by the full body weight. It is also possible to define a sub volume by picks or by numbers;

• Local SAR: SAR is given as a numerical value per volume element and becomes a space distribution function which is called SAR distribution. For this function, the mass mean value in arbitrary tissue volume is called local SAR.

For the local SAR calculation, one has to specify the mass in g over which the SAR should be averaged or choose point SAR. Typical values are 1 g or 10 g. For point SAR calculations, the averaging mass is denoted as zero. There are several methods to define an averaging volume which contains the desired averaging mass in CST. The IEEE C95.3 method uses a cube which is grown around the specific point P until it contains the required mass. The cubes are always aligned to the x yz coordinate system. However, if one face of the cube is completely in air, the volume is regarded as invalid and the SAR is taken from the point with the maximum SAR whose averaging volume contains the point. A modification of this method which is called CST C95.3 takes for invalid averaging volumes not the maximum but the closest valid SAR. The third method is IEEE/IEC 62704-1, and it is an averaging method according to the current dual logo standard draft IEEE/IEC 62704-1. Compared to the IEEE C95.3 method it uses an additional criterion that limits the air volume in valid averaging cubes to 10%.

In the present thesis, local SAR calculations are performed, and the maximum of the spatially averaged SAR for the given averaging mass are presented, for example, the peak 10 g average SAR.

Uncertainty in modelling

In [12], a brief history of modelling and the details of UQ in modelling are provided. Here, one for the most part borrows key elements from this description since very pertinent to the present endeavor.

Yet, models depend upon simplifying assumptions usually validated afterwards from experimental evidence. Moreover, they must be fed by values of input parameters, the estimation of which might be difficult or inaccurate, i.e., uncertain.

Fig. 1.1 Synoptic view of modelling

The framework of modelling is sketched in Fig. 1.1. It shows that there are modelling error, numerical error and uncertainty existing in modelling. Model error comes from lack of knowledge of the underlying true physics. It depends on how accurately a mathematical model describes the true system for a real-life situation, considering the fact that models are almost always only approximations to reality. Most mathematical models are too complicated to be tackled exactly. Some algorithms may be used to approximate the solution of the mathematical equations, which, however, introduces numerical errors, for example the Finite Element Method (FEM) [13] or Finite Difference Method (FDM) [14] may be used to approximate the solution of a partial differential equation. The definition of uncertainty is a potential deficiency in any phase or activity of the modelling process that is due to lack of knowledge. Uncertainty can then be classified into two categories [15]:

• "Type A" uncertainty or aleatory uncertainty which arises from stochastic or random behaviour in the system. It does not exist in deterministic computational simulations.

• "Type B" uncertainty or epistemic uncertainty which is the uncertainty that arises due to lack of knowledge. This type of uncertainty can be reduced by acquiring more knowledge. It can be further categorised into two groups. The first group Introduction of epistemic uncertainty is the uncertainty in how well the mathematical model represents the true behaviour of the real physical system. This uncertainty, which is known as model form uncertainty, is very difficult to determine but can be reduced by verifying the model against physical measurements. The second group of epistemic uncertainty is the uncertainty that arises due to a lack of precise input parameter data. For a system, the uncertainties in input parameters will induce the uncertainties in outputs. This type of uncertainty is often known as parameter uncertainty.

The framework of UQ is already presented in [12] as already said, and it is depicted in Fig. 1.2. Dealing with parameter uncertainty is the main goal of this thesis.

Fig. 1.2 The framework of uncertainty analysis

Three steps for uncertainty analysis have been identified:

• Step A consists in defining the model or sequence of models and associated criteria (e.g. failure criteria) to assess the physical system. For complex cases, it requires clear identification of the input and output of each sub-model. This step gathers all the ingredients used for a classical deterministic analysis of the physical system.

• Step B consists in quantifying sources of uncertainty, i.e., identifying those input parameters not well-known and modelling them within a probabilistic context. This provides a random vector of input parameters. In some cases, the description of the variability of parameters requires introduction of random processes.

• Step C consists in propagating the uncertainty in the input through the model, i.e., characterizing the random response appropriately vs. the assessment criteria of Step A. Numerous methods exist to do so.

• UQ methods provide information about the respective impact of the random input parameters on the response randomness. Such a hierarchization of the input parameters is known as sensitivity analysis. This type of analysis (Step C' in Fig. 1.2) might be the unique goal of the probabilistic study.

Step A (models of physical systems): a model of a real-world system is taken as a general function f : x → y = f (x). It can be explicitly written. More often, it is a black box function like a computer program that takes input values and yields a result. Then, the model f is known only through point-wise evaluations y i = f (x i ) for each input vector x i for which the program is run. The input parameters can be geometrical ones (dimensions, length of structural members, cross-sections, etc.), material ones (permittivity, permeability and conductivity), and parameters describing the boundary conditions. The output of a model can be many things obviously. SAR values, e.g., peak 10g average and peak 1g average are the main focus here, but fields and currents within a discretized version of the model and responses of probes under many forms can be of interest as well. Also, even though the model is sufficiently accurate (equations relevant to describe the physical phenomena, approximations in the computational scheme acceptable), the input parameters are not perfectly known, and uncertainty in input parameters exists.

Step B (quantification of sources of uncertainty): A real system never fully matches its initial design. True dimensions do not correspond exactly due to imperfections in manufacturing. Material properties differ slightly from the codified properties of the class of material which it is supposed to be made of. Loading (i.e., boundary conditions) is idealized, only roughly representing the complexity of the one the real system. As a consequence, modelling input parameters in a probabilistic context makes sense. The end product of Step B is the description of the random vector of input parameters X in terms of joint Probability Density Function (PDF) W X (x).

Step C (uncertainty propagation): considering a model x → y = f (x), uncertainty propagation methods are required to study the probabilistic content of the output. Three categories of methods appear. Second moment ones compute the mean value and variance of the model response (higher-order moment methods may be envisaged). Structural reliability ones essentially investigate the tails of the response PDF by getting the probability of exceeding a threshold (probability of failure). Spectral ones (stochastic Introduction finite element methods in computational mechanics) represent the complete randomness of the response in intrinsic fashion via functional analysis. Problems pertaining to second moment or structural reliability analysis can be solved by post-processing their basic output, i.e., expansion coefficients on a suitable basis of functions.

Step C' (sensitivity analysis): the modelling of complex systems usually requires many input parameters. Defining a proper probabilistic model for each one is a burden. However, in real-world, only a limited number of input parameters tends to influence the response randomness. Sensitivity analysis aims at selecting the important parameters, usually via quantitative importance measures or sensitivity indices. Most methods for uncertainty propagation provide sensitivity measures as a by-product of their principal output.

Probability approach for uncertain quantification

The general framework of uncertainty analysis as given in Fig. 1.2, shows that the analyst has to build a probabilistic description of the data and more precisely, of the input parameters of the model. The input parameters of a model can be modelled as random variables or random vectors.

Random variables

A real random variable is a mapping X : Ω -→ D X ⊂ . When D X is a discrete set, the random variable is said discrete, otherwise it is said continuous. For continuous random variables, they have continuous range of values within a specific interval. In this case, an exact value of the random variable cannot be found. This is the reason why intervals are used to ascertain the likelihood that a random variable is between two known values. The physical parameters of a numerical simulation model are continuous random variables. The discrete random variables arise from discrete values of a set. A random variable is completely defined by its Cumulative Distribution Function (CDF) denoted by F X (x):

F X (x) = (X ≤ x) (1.6)
where is the probability. For a discrete random variable, the domain of definition may be represented as D X = {x (i) , i ∈ }. The probability mass function is defined in this case as follows:

p i = (X = x (i) )
(1.7)

1.2 Probability approach for uncertain quantification 9

The CDF is:

F X (x) = i∈ p i S {x≥x (i) } (x) (1.8)
where S {x≥x (i) } (x) is the indicator function of the set x ∈ : x ≥ a defined by:

S x≥a (x) = 1 i f x ≥ a 0 other wise (1.9)
For a continuous random variable, the probability density function (PDF) is defined as:

W X (x) = lim m→0,m>0 (x ≤ X ≤ x + h)/h (1.10)
and

W X (x) = d F X (x) d x (1.11)
The mathematical expectation is denoted by E[.], and the mean is defined as the first moment about the origin which acts as a representative of the central value for the distribution. A mathematical expression for the mean denoted by µ X of a continuous random variable with a probability density function PDF W X (x) is:

µ X = E[X ] = D X xW X (x)d x (1.12)
The expectation of a function h(X ) is defined as:

E[h(x)] = D X h(x)W X (x)d x (1.13)
The variance denoted by σ 2 X , or its square root, the standard deviation, is a measure available for the second central moment of data and is given by:

σ 2 X = E[(X -µ X ) 2 ] = D X (x -µ X ) 2 W X (x)d x (1.14)
The third central moment about the mean describes the symmetry or skew of the distribution. The skew for a continuous random variable is defined as follows:

γ 1 = E X -µ X σ X 3 = D X x -µ X σ X 3 W X (x)d x (1.15)

Introduction

The fourth central moment about the mean describes the relative peakedness or flatness of a probability density function relative to a normal distribution. This kurtosis is a non-dimensional quantity related to the fourth central moment. The kurtosis for a continuous random variable is defined by the expression:

γ 2 = E X -µ X σ X 4 -3 = D X x -µ X σ X 4 W X (x)d x -3
The covariance of two random variables X and Y is:

C ov[X , Y ] = E[(X -µ X )(Y -µ Y )] (1.16)
In statistics, dependence is any statistical relationship between two random variables or two sets of data. Correlation refers to any of a broad class of statistical relationships involving dependence. There are several correlation coefficients, often denoted ρ X Y , measuring the degree of correlation. The most common one is the Pearson correlation coefficient [START_REF] Menegaz | A new smallest sigma set for the Unscented Transform and its applications on SLAM[END_REF], which is sensitive only to a linear relationship between two variables. Other correlation coefficients have been developed to be more robust than the Pearson correlation, and more sensitive to non-linear relationships. The "Pearson's correlation coefficient" is commonly called simply "the correlation coefficient", and it is:

ρ X Y = cor r(X , Y ) = C ov(X , Y ) σ X σ Y = E[(X -µ X )(Y -µ Y )] σ X σ Y (1.17)
The Pearson correlation is +1 in the case of a perfect direct linear relationship (correlation), -1 in the case of a perfect decreasing linear relationship, and the value between -1 and 1 in all other cases, indicating the degree of linear dependence between the variables. As it approaches zero there is less of a relationship (closer to uncorrelated).

The closer the coefficient is to either -1 or 1, the stronger the correlation between the variables. If the variables are independent, Pearson's correlation coefficient is 0, but the converse is not true because the correlation coefficient only detects linear dependencies between two variables [START_REF] Menegaz | A new smallest sigma set for the Unscented Transform and its applications on SLAM[END_REF].

It is always possible to remove the correlation between random variables with a linear transformation, even if the relationship between the variables is non-linear [START_REF] Menegaz | A new smallest sigma set for the Unscented Transform and its applications on SLAM[END_REF]. A presentation of this result for population distributions is given by Cox and Hinkley [17].

Gaussian random variables are of utmost importance in probability theory. They are denoted by X ∼ (µ X , σ X ) where µ X is the mean value and σ X is the standard deviation. The PDF W X (x) is defined by:

W X (x) = 1 2πσ 2 X e - (x-µ X ) 2 2σ 2 X (1.18)
The standard normal CDF F X (x) is:

F X (x) = x -∞ 1 2πσ 2 X e - (x-µ X ) 2 2σ 2 X d x (1.19)
When µ X = 0, σ X = 1,the standard normal PDF W X (x) is defined by:

W X (x) = 1 2π e -x 2 /2 (1.20)
The standard normal CDF F X (x) is:

F X (x) = x -∞ 1 2π e -x 2 /2 d x (1.21)

Random vectors

A real random vector X is a mapping X : Ω -→ D X ⊂ N , where N is the dimension or size of the vector and N ≥ 2. The random vector X ≡ {X 1 , . . . , X N } T has components, and these components are considered as random variables. (.) T denotes the transposition. The joint PDF is denoted by W X (x), and the marginal distribution of a given component X i is obtained by integrating the joint PDF over all the remaining components:

W X i (x i ) = D X ∼i W X (x)d x ∼i (1.22)
where

d x ∼i = d x 1 . . . d x i . . . d x N , and
D X ∼i is the subset of D X defined by {x ⊂ D X }.
Similarly, the joint distribution of two components (X i , X j ) is given by:

W X i ,X j (x i , x j ) = D X ∼(i, j) W X (x)d x ∼i j (1.23)
The expectation of a random vector is the vector containing the expectation of each component:

µ X = {µ X 1 , . . . , µ X N } T (1.24)

Introduction

The covariance matrix of X is a square symmetric matrix C of size N , given by:

C = ΛRΛ (1.25)
where

C i, j = C ov[X i , X j ] (1.26) R i, j = ρ X i ,X j (1.27)
and the diagonal matrix Λ = Diag(σ X 1 , . . . , σ X N ) contains the standard deviation of each component of X.

Gaussian random vectors denoted by X ∼ (µ X , C), are completely defined by their mean value vector µ X and covariance matrix C through the following joint PDF:

W X (x) = (2π) -N 2 (detC) -1 2 ex p - 1 2 (x -µ X ) T C -1 (x -µ X ) (1.28)
Standard normal random vectors of size N are defined by X ∼ (0, I N ) where I N is the identity matrix. The PDF is:

W X (x) = (2π) -N 2 ex p - 1 2 (x 2 1 + • • • + x 2 N ) (1.29)

Overview of the methods related to uncertainty quantification

Uncertainty analysis has been studied in many areas such as mechanical engineering [12], climate modelling [18], Computational Fluid Dynamics (CFD) [19], etc. In the field of Computational Electromagnetics (CEM), an IEEE standard: IEEE P1597.1 [20] has been approved for the "Validation of Computational Electromagnetics Computer Modeling and Simulations". This standard outlines the necessary steps for validating CEM models and simulations, and much work carried out by more and more researchers to provide useful contributions for future editions of the IEEE standard. Future editions of this standard should require that the results of all CEM simulations are accompanied by an approximate estimate of the uncertainty in the results [21]. In 2006, the publication of Chauvière [START_REF] Chauvière | Computational modelling of uncertainty in time-domain electromagnetics[END_REF] presented the work of implementing the Polynomial Chaos (PC) method into a higher-order discontinuous Galerkin solution of Maxwell's equations [START_REF] Chauvière | Computational modelling of uncertainty in time-domain electromagnetics[END_REF]. In 2008, the use of a Direct Solution Technique (DST) [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF] to quantify uncertainty 1.3 Overview of the methods related to uncertainty quantification 13

was presented by Ajayi . In his thesis, he went on estimating the uncertainty in the first resonant frequency of simple electromagnetic problems using the Method of Moments (MoM) [12, 21, 23][24-26]. The use of Unscented Transforms (UT) [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF][START_REF] Julier | Consistent debiased method for converting between polar and Cartesian coordinate systems[END_REF][START_REF] Menegaz | A new smallest sigma set for the Unscented Transform and its applications on SLAM[END_REF][START_REF] De Menezes | Efficient computation of stochastic electromagnetic problems using unscented transforms[END_REF] to efficiently estimate uncertainties in CEM simulation has also been proposed in [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF][START_REF] De Menezes | Efficient computation of stochastic electromagnetic problems using unscented transforms[END_REF].

In 2009, in his thesis [21], Edwards presented how to implement PC expansion into the FDTD code for UQ in CEM simulation. The author considered three uncorrelated uncertain inputs in the simulation.

The methods for propagating the uncertainty in the input parameters through a model may be classified according to the type of information that is investigated with respect to the output of the model [12]. In this reference, the author proposed the following classification of methods for UQ:

• When the mean value and standard deviation of the response are of interest, second moment analysis methods are to be used. Methods such as Monte Carlo method, the perturbation method (MoM), the weighted integral method or the quadrature method enter into this category. They provide a first estimate of the response variability, which is limited to the first two statistical moments of the response.

• When the tail of the response PDF is of interest (i.e., the low quantiles), methods of structural reliability analysis such as importance sampling or directional simulation may be used in this case.

• When the whole PDF is of interest, methods of approximation thereof have to be considered. Monte Carlo Simulation (MCS) is the basic approach to solve the problem.

Since only the mean value and standard deviation of the response are of interest in this thesis, the UQ methods used herein are both second moment analysis methods. Below, one important and popular second moment analysis method (MoM), and the quadrature method which is used in some UQ methods are presented.

Method of moments

The Method of Moments (MCM) is also known as perturbation method which is a nonsampling method based on the Taylor series expansion of the model response around the mean value of the input parameters [12]. Typically, a second-order expansion is employed because the resulting system of equations becomes extremely complicated beyond the second order. This method is called MoM throughout this thesis.

Introduction

In [21], the author indicates that MoM is the internationally accepted method outlined in the Guide to the Expression of Uncertainty in Measurement (GUM) [START_REF]Joint Committee for Guides in Metrology, Evaluation of measurement data -Guide to the expression of uncertainty in measurement[END_REF] for the propagation of uncertainties through a model. A similar method is outlined in the United Kingdom Accreditation Service (UKAS) [15] for the estimation of the uncertainty in practical Electromagnetic Compatibility (EMC) measurements. In [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF], the MoM is applied to a two-dimensional Transmission Line Matrix (TLM) [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF] method, and is shown to be accurate for small parameter variations, whilst being computationally more efficient than the MCM. The MoM calculates the mean and uncertainty of the responses of a model by constructing truncated Taylor series expansions around the expected value of the input parameters [START_REF] Walters | Uncertainty analysis for fluid dynamics with applications[END_REF].

The idea of the MoM is [21]: considering f be a function of only one variable denoted by x, the output is y = f (x). f can be expanded around the mean input value µ x to give

f (x) = f (µ x ) + d f d x x=µ x (x -µ x ) + d 2 f d x 2 x=µ x (x -µ x ) 2 2 + . . . (1.30)
Taking the expected value of both sides of equation (1.30) yields

E[ y] = E[ f (x)] = f (µ x ) + 1 2 d 2 f d x 2 x=µ x σ 2 x + . . . (1.31)
where σ 2 x is the variance in x. The first order estimate of the mean µ y1st of the output y is given by

µ y1st = f (µ x ) (1.32)
and it is obtained by truncating the initial Taylor series to first order.

The second-order estimate of the mean µ y2nd is obtained by truncating the Taylor series to second order, and it is given by

µ y2nd = f (µ x ) + 1 2 d 2 f d x 2 x=µ x σ 2 x (1.33)

Overview of the methods related to uncertainty quantification 15

The variance in the output y is calculated as follows:

E[ y 2 ] = E[ f 2 (x)] = f 2 (µ x ) + 1 2 d 2 ( f 2 ) d x 2 x=µ x σ 2 x = f 2 (µ x ) + d f d x 2 + f d 2 f d x 2 x=µ x σ 2 x (1.34)
and then the first order estimate of the uncertainty in y denoted by σ y1st , is given by

σ 2 y1st = E[ y 2 ] -E[ y] 2 = f 2 (µ x ) + d f d x 2 + f (µ x ) d 2 f d x 2 x=µ x σ 2 x -f (µ x ) + 1 2 d 2 f d x 2 x=µ x σ 2 x 2 = σ 2 x d f d x x=µ x 2 (1.35)
where terms involving σ 4

x are ignored. The second order estimate of the variance is:

σ 2 y2nd = σ 2 x d f d x x=µ x 2 - 1 2 d 2 f d x 2 x=µ x σ 2 x 2 (1.36)
These estimates of the uncertainty are formed using truncated Taylor series expansions, and they will therefore only be accurate if the input uncertainty σ x is small [21].

A brief description of the application of the method to FDTD simulations is given in [21]. Considering the electric field E as the output of an FDTD simulation. The output can be given as a function of the frequency f and N uncertain input parameters x 1 , . . . , x N :

E = E[ f , x 1 , . . . , x N ] (1.37)
To calculate the variance in the output due to the variance in the input parameters, the sensitivity derivative of the output with respect to each individual input parameter must first be found. The sensitivity of the output to the parameter △x i is approximated using a forward difference approximation to give

d E d x i f ,µ 1 ,...,µ N = E( f , µ 1 , . . . , µ i + △x i , . . . , µ N ) -E( f , µ 1 , . . . , µ i , . . . , µ N ) △x i (1.38)
where {µ 1 , . . . , µ N } is the mean of {x 1 , . . . , x N } and △x i is a small change in the ith input parameter x i . For a FDTD simulation this amounts to performing one simulation with the input parameters taking on their mean values to obtain E( f , µ 1 , . . . , µ N ), and one simulation with the ith parameter value perturbed by a small value △x i . To obtain all N sensitivity derivatives, which are required to calculate the variance in E, N + 1 simulations must be performed. The computational runtime should however be significantly less than the one required by the MCM [21]. The limitation of the MoM is that the magnitude of the uncertainty at both input parameters and outputs cannot be too large, otherwise the method does not give accurate predictions.

Quadrature method

The quadrature method [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] is a numerical integration technique for an integral of the form:

D X f (x)W X (x)d x (1.39)
and the quadrature rule can be stated as:

D X f (x)W X (x)d x = N i=1 w i f (x i ) (1.40)
where w i are constants which depend only on the weighting function W X (x) and x i are constants in the region of integration. Using the quadrature method, the mean and variance of the response given by:

µ y = E[ y] = D X f (x)W X (x)d x (1.41) σ 2 y = E[( y -µ y ) 2 ] = D X [ f (x) -µ y ] 2 W X (x)d x (1.42)
can be evaluated. Assume now that X is a random vector with independent components, i.e., its PDF is:

W X (x) = W X 1 (x 1 ) . . . W X N (x N ) (1.43)
1.4 Aims and objectives 17

D X f (x)W X (x)d x = m i 1 j 1 =1
...

m i N j N =1 f (x i 1 j 1 , ..., x i N j N )(w i 1 j 1 ⊗ ...⊗w i N j N ) (1.44)
where multi-index {i 1 , i 2 , ..., i N } T ∈Z + N , and ⊗ denotes tensor product. It is obvious that the above product needs N n=1 m i n function evaluations. The number of collocation points in the tensor grid grows exponentially with the number of input random variables. Using full tensor product quadrature method for estimation faces the so-called curse of dimensionality.

Aims and objectives

The main aim of this thesis is to quantify the parameter uncertainty in SAR calculation of a mobile phone. MCM is a traditional UQ method but it is not efficient in SAR calculation because it converges slowly, resulting in high computational cost. A typical SAR calculation of a mobile phone placed against the Specific Anthropomorphic Mannequin (SAM) head phantom [START_REF] Beard | Comparisons of computed mobile phone induced SAR in the SAM phantom to that in anatomically correct models of the humanhead[END_REF] requires about half an hour using a Graphical Processing Unit (GPU) to several hours using a standard workstation.

Alternative UQ methods must be investigated for this challenging work. The methods used for UQ can be either intrusive or non-intrusive. The application of intrusive UQ methods requires modifying the numerical code in the simulation solver. On the other hand, using non-intrusive UQ methods, the simulation process is simply considered as a black-box without changing the code of the simulation solver. However, intrusive UQ methods can provide the statistical information of the outputs by one simulation instead of running the code several times.

In the present work, three non-intrusive UQ methods are investigated for UQ in SAR calculation of a mobile phone. These methods are: UT, Stochastic Collocation (SC) [START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF]35] and Non-Intrusive Polynomial Chaos (NIPC) [START_REF] Walters | Uncertainty analysis for fluid dynamics with applications[END_REF].

UQ methods share a common disadvantage due to the fact it is difficult to use them with multiple uncertain inputs in the calculation. The formulation of multidimensional UQ method is complex, and the number of simulations required for UQ increases with the number of input parameters. However, in SAR calculation there are usually several random variables due to the lack of precise knowledge of the material or geometry of the mobile phone being modelled. To solve this problem, in this thesis, combined uncertainties [START_REF]Joint Committee for Guides in Metrology, Evaluation of measurement data -Guide to the expression of uncertainty in measurement[END_REF] in SAR calculation is proposed.

A kind of variance-based Sensitivity Analysis (SA) method: PC based Sobol' indices [12,[START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF][START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] is introduced for SA in SAR calculation. The characteristic of PC based Introduction Sobol' indices is that the sensitivity indices are calculated by PC expansion. When the uncertainty in SAR simulation is quantified by NIPC, the sensitivity indices can be obtained by a further calculation following the UQ calculations. UQ and SA can be implemented simultaneously.

Outline of thesis

• Chapter 2 first introduces the traditional UQ method: MCM. Three sampling methods applied in MCM are presented, and they are random sampling, stratified sampling and Latin Hypercube Sampling (LHS). Then the UT method which is a similar method to MCM is described in detail. The idea of UT, the formulation of UT with one or two random variables and the form of high-order UT are given. Finally, the results of UT and MCM are compared for three examples using a simple mobile phone model.

• Chapter 3 presents two similar methods: SC expansion and PC expansion. Using UT, SC and NIPC for UQ in SAR calculation is investigated in case of one or two random variables in the calculation, and the results of these three methods are given and analysed. One of variance-based SA methods: PC based Sobol' indices, is described and used in the calculation. The idea of combined uncertainties is presented, and the method is validated in SAR calculation. A conclusion follows.

• Chapter 4 first introduces the stochastic system formulation which is the foundation of intrusive PC expansion. Then the one-dimensional FDTD (1D FDTD) method and gPC expansion in 1D FDTD are presented. A conclusion about intrusive PC follows.

• Chapter 5 provides the major conclusions, and future continuation of this work is also considered.

Chapter 2

Uncertainty quantification using the Unscented Transformation (UT)

Introduction

MCM is the most recognized method to estimate the uncertainty in various CEM simulations. The main drawback of MCM is that it converges slowly. It is difficult to implement MCM in SAR calculations due to this drawback, because a typical SAR calculation of a mobile phone placed against the SAM head phantom requires about half an hour or even several hours depending on the computing machine used for the simulation. It is not practical to perform several thousands SAR calculations, and some alternative methods should be proposed for UQ in SAR calculation. In this chapter, one of the computationally efficient UQ methods is proposed for UQ in SAR calculation. This method is referred to as UT. In essence, UT is rather similar to MCM since both are used to find an approximation of the probability distribution function rather than the non-linear function which is thus seen as a black box [START_REF] Julier | Consistent debiased method for converting between polar and Cartesian coordinate systems[END_REF]. UT is a non-intrusive UQ method which is proposed in [START_REF] Julier | Consistent debiased method for converting between polar and Cartesian coordinate systems[END_REF] and then used for UQ in the field of CEM [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF][START_REF] De Menezes | Efficient computation of stochastic electromagnetic problems using unscented transforms[END_REF].

The idea of UT is presented in detail in this chapter, and the results of applying UT and MCM for UQ in calculations are compared and analysed.

The Monte Carlo method

The name of MCM refers to the Monte Carlo casino in Monaco and was given to the method because of the use of probability and repetition [21]. The MCM works like this: the "numerical" experiment is modelled as an unknown transformation f (x) of the inputs x, which have a known PDF W x (x). It is obvious that a number of sets of random samples from the input space should be taken, and the random samples x 1 , . . . , x N of x must be selected as successive inputs sets in order to obtain the desired statistics of y.

The number of samples of random variable N should not be too large because of the running time of the code. Let x = {x 1 , . . . , x N } be a random vector, and the output is

y = f (x)
where f is a function of random variable x. The mean and variance of the output obtained by MCM are

µ y = 1 N N i=1 f (x i ) (2.1)
σ 2 y = 1 N -1 N i=1 ( f (x i ) -µ y ) 2 (2.2)
The sampling method which is applied in MCM influences the convergence speed of MCM [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. There are three main sampling methods: random sampling, stratified sampling and LHS. LHS is the preferred sampling method because compared with other two sampling methods, LHS lets the MCM converge more quickly. These three sampling methods are presented in this chapter.

Random sampling and stratified sampling

Choosing the input values x 1 , . . . , x N as random samples from the PDF W x (x), it is random sampling. Each individual is chosen randomly and entirely by chance. New sample points of random sampling are generated without taking into account the previously generated sample points. One does not necessarily need to know beforehand how many sample points are needed.

The stratified sampling is that all areas of the sample space of x are represented by input values. Let the sample space S of x be partitioned into I disjoint strata S i . Obtain a random sample x i j , j = 1, . . . , n i from S i . The n i sum to N , and if I = 1, we have random sampling over the entire sample space. Comparing with random sampling, the advantage of stratified sampling is that if population density varies greatly within a region, stratified sampling will ensure that estimates can be made with equal accuracy in different parts of the region, and that comparisons of sub-regions can be made with equal statistical power. The disadvantage of stratified sampling is it is not useful when the population cannot be exhaustively partitioned into disjoint subgroups. It would be a misapplication of the technique to make subgroups' sample sizes proportional to the amount of data available from the subgroups, rather than scaling sample sizes to subgroup sizes. Data to the subgroups' sizes within the total population. For an efficient way to partition sampling resources among groups that vary in their means, their variances, and their costs, see "optimum allocation".

Latin hypercube sampling (LHS)

LHS is proposed by McKay et al. [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF], and then it was further elaborated in 1977 [START_REF] Eglajs | New approach to the design of multifactor experiments[END_REF] and in 1981 [START_REF] Iman | An approach to sensitivity analysis of computer models, Part 1. Introduction, input variable selection and preliminary variable assessment[END_REF]. In the context of statistical sampling, a square grid containing sample positions is a Latin square if (and only if) there is only one sample in each row and each column. A Latin hypercube is the generalisation of this concept to an arbitrary number of dimensions, whereby each sample is the only one in each axis-aligned hyperplane containing it. Let x 1 , . . . , x K are K uncertain inputs. The PDF of the kth uncertain input x k is split into N segments of equal probability. One random sample is taken from each segment, producing N samples for parameter k that span the whole of the PDF. This process is repeated for all K parameters. By then randomly combining the N samples from the K parameters, N sets of input parameter values are produced. The maximum number of combinations for a Latin hypercube of N divisions and K variables (i.e., dimensions) can be computed with the following formula:

N -1 i=0 (N -i) K-1 = (N !) K-1 (2.3)
For example, a Latin hypercube of N = 4 divisions with K = 2 variables will have 24 possible combinations. In LHS, the number of sample points and segments must be decided firstly. One advantage of the LHS is when the output y = f (x) is dominated by only a few of the components of x, and it ensures that each of those components is represented in a fully stratified manner, no matter which components might turn out to be important. Using LHS with the MCM has been shown to increase the rate of convergence compared to using random sampling or stratified sampling [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF].

Unscented transformation (UT)

UT is proposed by Julier and Uhlmann in 1997 [START_REF] Julier | Consistent debiased method for converting between polar and Cartesian coordinate systems[END_REF] as a method to propagate mean and covariance information through non-linear transformations. The authors' motivation for proposing UT is that applying the Extended Kalman Filter (EKF) [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF] to non-linear systems is difficult to implement, difficult to tune, and only reliable for systems that are Uncertainty quantification using the Unscented Transformation (UT) almost linear on the time scale of the updates, so an alternative must be proposed as a practical estimator to non-linear systems. Then the use of UT to efficiently estimate uncertainties has been proposed in [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF][START_REF] De Menezes | Efficient computation of stochastic electromagnetic problems using unscented transforms[END_REF].

The basic idea of UT

In [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF], the author gives detailed description about UT. Considering a non-linear mapping f of a random variable x with mean µ x and standard deviation σ x , the output y of the non-linear transformation is simply expressed as:

y = f (x) (2.4)
and the objective is to obtain the statistics of y, such as the mean µ y and variance σ 2 y , from which the standard deviation is derived. For this objective, let x be a random variable with zero mean and the same standard deviation as x, that is x = 0 ± σ x . Respective to this, x is expressed as:

x = µ x + x (2.5)
The Taylor's series approximation for a function of x is:

f (x) = f (µ x + x) = f (µ x ) + ( d f d x | x=µ x )x + 1 2! ( d 2 f d x 2 | x=µ x )x 2 + 1 3! ( d 3 f d x 3 | x=µ x )x 3 + . . . (2.6)
and equation (2.6) can be simplified into the following form:

f (µ x + x) = a 0 + a 1 x + a 2 x2 + • • • + a n x n = a 0 + h(x) (2.7)
where

a 0 = f (µ x ), a n = 1 n! ( d n f d x n | x=µ x )
x n and n is the order of the polynomial. The expected value of f (x) is:

µ y = E[ f (x)] = E[ f (µ x + x)] = a 0 + E[h(x)] = f (µ x ) + µ h (2.8)
The variance of σ 2 y is:

σ 2 y = E[( f (x) -µ y ) 2 ] = E[h 2 (x)] -µ h 2 (2.9)
The sigma points S i used in the UT calculations are defined at specific values and designated with weights w i . As those are discrete values of stochastic variables, the expected value of f (x) in terms of these parameters is

µ y ≈ w 0 f (µ x + S 0 ) + m-1 i=1 w i f (µ x + S i ) (2.10)
where m is the number of sigma points, and S 0 = 0. Substituting (2.7) into (2.10), we obtain the following equation (keeping " = " signal for simplest, though from equation (2.10) approximation are used.):

µ y = [w 0 + m-1 i=1 w i ] f (µ x ) + m-1 i=1 w i h(S i ) (2.11)
Considering equations (2.8) and (2.11), the following relations are obtained:

[w 0 + m-1 i=1 w i ] = 1 (2.12) m-1 i=1 w i h(S i ) = µ h (2.13)
Equation (2.12) shows that the weights w i in UT are deterministic and must sum up to one. Similarly, the variance σ 2 y in terms of weights and sigma points is

σ 2 y = w 0 [ f (µ x ) -µ y ] 2 + m-1 i=1 w i [ f (µ x + S i ) -µ y ] 2 (2.14)
Substituting (2.8) into (2.14), we obtain the following equation:

σ 2 y = w 0 µ h 2 + m-1 i=1 w i [h(S i ) -µ h ] 2 (2.15)
Expanding the above equation:

σ 2 y = [w 0 + m-1 i=1 w i ]µ h 2 + m-1 i=1 w i [h 2 (S i ) -2µ h h(S i )] (2.16)
Uncertainty quantification using the Unscented Transformation (UT) Substituting (2.12) and (2.13) into (2.16):

σ 2 y = m-1 i=1 [h 2 (S i )] -µ h 2 (2.17)
Comparing (2.9) and the equation (2.17) gives the relation:

E[h 2 (x)] = m-1 i=1 w i h 2 (S i ) (2.18)
In general, the sigma points and weights in UT are related to the moments of a random variable x by:

E[x k ] = m-1 i=1 w i S k i (2.19)
where k is the order of the expansion. For example, a second-order UT (n = 2) results in a non-linear system with four variables (w 1 , w 2 , S 1 , S 2 ). Since it requires four equations to obtain these four variables, setting k = 4, then the first four moments are obtained.

The first four moments of any distribution provide information on the mean, variance, skew and kurtosis. From these moments the following sets of equations are obtained:

E[x] = m-1 i=1 w i S i = 0 (2.20) E[x 2 ] = m-1 i=1 w i S 2 i = σ 2 x (2.21) E[x 3 ] = m-1 i=1 w i S 3 i = γ 1 σ 3 x (2.22) E[x 4 ] = m-1 i=1 w i S 4 i = (γ 2 + 3)σ 4 x (2.23)

Unscented transformation (UT)

25

where γ 1 and γ 2 are the skew and kurtosis of the input distribution. Solving equations (2.20)-(2.23), the sigma points S i and weights w i can be obtained as follows:

S 1 = 1 2 γ 1 -4(γ 2 + 3) -3γ 2 1 σ x S 2 = 1 2 γ 1 + 4(γ 2 + 3) -3γ 2 1 σ x w 1 = -2 4(γ 2 + 3) -3γ 2 1 γ 1 -4(γ 2 + 3) -3γ 2 1 w 2 = 2 4(γ 2 + 3) -3γ 2 1 γ 1 + 4(γ 2 + 3) -3γ 2 1 (2.24)
In the case of a Gaussian distribution, the skew γ 1 = 0 and kurtosis γ 2 = 0 result into the sigma points 3σ x and -3σ x and weights w 0 = 2 3 and w 1 = w 2 = 1 6 for second-order UT. If the same order is used with the uniform distribution, the results are

S 1 = -3 1 5 σ x , S 2 = 3 1 5 σ x , w 1 = 5
78 and w 2 = 5 78 with the γ 1 = 0 and γ 2 = -6 5 . The mean µ y and variance σ 2 y of UT are

µ y = m-1 i=0 w i f (µ x + S i ) (2.25) σ 2 y = m-1 i=0 w i ( f (µ x + S i ) -µ y ) 2 (2.26)
The weights w i and sigma points S i can also be extended to find the skew γ 1 y and kurtosis γ 2 y of the output y:

γ 1 y = m-1 i=0 w i f (µ x + S i ) -µ y 3 σ 3 y (2.27) γ 2 y = m-1 i=0 w i f (µ x + S i ) -µ y 4 σ 4 y -3 (2.28)

The formulation of two-dimensional UT

When there are two random variables x 1 and x 2 with x1 and x2 respectively in the problem, the Taylor polynomial representation is still suitable and it reads in terms of Uncertainty quantification using the Unscented Transformation (UT) regression coefficients as [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF]:

h(x 1 , x2 ) = a 1 x1 + a 2 x2 + a 3 x1 x2 + a 4 x2 1 + a 5 x2 2 + . . . (2.29)
Equation (2.19) can be modified, and the relationship between sigma points and moments is:

E[x k 1 x l 2 ] = m-1 i=1 w i (S k 1i S l 2i ) (2.30)
In the case of the truncated second-order (n = 2) polynomial which is given by

h(x 1 , x2 ) = a 1 x1 + a 2 x2 + a 3 x1 x2 + a 4 x2 1 + a 5 x2 2 (2.31)
the combination of k and l, 1 ≤ k + l ≤ 4, leads to the system:

m-1 i=1 w i S 1i = E[x 1 ] m-1 i=1 w i S 2i = E[x 2 ] m-1 i=1 w i S 2 1i = E[x 2 1 ] m-1 i=1 w i S 2 2i = E[x 2 2 ] m-1 i=1 w i S 3 1i = E[x 3 1 ] m-1 i=1 w i S 3 2i = E[x 3 2 ] m-1 i=1 w i S 4 1i = E[x 4 1 ] m-1 i=1 w i S 4 2i = E[x 4 2 ] m-1 i=1 w i S 1i S 2i = E[x 1 x2 ] m-1 i=1 w i S 1i S 2 2i = E[x 1 x2 2 ] m-1 i=1 w i S 2 1i S 2i = E[x 2 1 x2 ] m-1 i=1 w i S 2 1i S 2 2i = E[x 2 1 x2 2 ] m-1 i=1 w i S 3 1i S 2i = E[x 3 1 x2 ] m-1 i=1 w i S 1i S 3 2i = E[x 1 x3 2 ]
(2.32)

There are 14 equations in (2.32) and each sigma point contributes with three variables w i , S 1i and S 2i . Therefore, five sigma points is the minimum amount that satisfies all conditions described in (2.32). The system has now 14 equations and 15 variables.

There are several possible solutions. In the case of zero mean independent normal random variables one solution is [START_REF] De Menezes | Efficient computation of stochastic electromagnetic problems using unscented transforms[END_REF]:

w 1 = w 2 = w 3 = w 4 = w 5 = 1 10 w 0 = 1 -(w 1 + w 2 + w 3 + w 4 + w 5 ) S 11 = 2cos(0 2π 5 )σ x 1 S 21 = 2sin(0 2π 5 )σ x 2 S 12 = 2cos( 2π
5 )σ x 1 S 22 = 2sin( 2π
5 )σ x 2 S 13 = 2cos(2 2π
5 )σ x 1 S 23 = 2sin(2 2π
5 )σ x 2 S 14 = 2cos(3 2π
5 )σ x 1 S 24 = 2sin(3 2π
5 )σ x 2 S 15 = 2cos(4 2π
5 )σ x 1 S 25 = 2sin(4 2π
5 )σ x 2 (2.33)
where σ x 1 and σ x 2 are standard deviations of x1 and x2 . The equations mentioned above are used for two independent random variables. If the two random variables are correlated such as they are two correlated normally distributed random variables, the covariance matrix transformation must be used as follows [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF]:

    S 1 σ x 1 S 2 σ x 2     =    1+ρ 12 2 1-ρ 12 2 1+ρ 12 2 - 1-ρ 12 2       S 1 σ x 1 S 2 σ x 2    (2.34)
Then one has the following new sigma points:

S 1 = σ x 1 1 + ρ 12 2 S 1 σ x 1 + 1 -ρ 12 2 S 2 σ x 2 (2.35) S 2 = σ x 2 1 + ρ 12 2 S 1 σ x 1 - 1 -ρ 12 2 S 2 σ x 2 (2.36)
where S 1 and S 2 are the sigma points for the two correlated random variables. ρ 12 is the correlation coefficient between x 1 and x 2 . The solutions of the equations are as Uncertainty quantification using the Unscented Transformation (UT)

follows [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF]:

w 1 = w 2 = w 3 = w 4 = w 5 = 1 10 w 0 = 1 -(w 1 + w 2 + w 3 + w 4 + w 5 ) S 11 = 2 A 1 cos(0 2π 5 )σ x 1 + A 2 sin(0 2π 5 )σ x 1 S 12 = 2 A 1 cos( 2π 5 )σ x 1 + A 2 sin( 2π 5 )σ x 1 S 13 = 2 A 1 cos(2 2π 5 )σ x 1 + A 2 sin(2 2π 5 )σ x 1 S 14 = 2 A 1 cos(3 2π 5 )σ x 1 + A 2 sin(3 2π 5 )σ x 1 S 15 = 2 A 1 cos(4 2π 5 )σ x 1 + A 2 sin(4 2π 5 )σ x 1 S 21 = 2 A 1 cos(0 2π 5 )σ x 2 -A 2 sin(0 2π 5 )σ x 2 S 22 = 2 A 1 cos( 2π 5 )σ x 2 -A 2 sin( 2π 5 )σ x 2 S 23 = 2 A 1 cos(2 2π 5 )σ x 2 -A 2 sin(2 2π 5 )σ x 2 S 24 = 2 A 1 cos(3 2π 5 )σ x 2 -A 2 sin(3 2π 5 )σ x 2 S 25 = 2 A 1 cos(4 2π 5 )σ x 2 -A 2 sin(4 2π 5 )σ x 2 (2.

37)

where

A 1 = 1 + ρ 12 2 A 2 = 1 -ρ 12 2
The mean and the variance of output for two-dimensional UT are as follows:

µ y = n i=0 n j=0 w i w j f (µ x 1 + S i , µ x 2 + S j ) (2.38) σ 2 y = n i=0 n j=0 w i w j ( f (µ x 1 + S i , µ x 2 + S j ) -µ y ) 2 (2.39) 2.
3 Unscented transformation (UT)

29

where x 1 and x 2 are two random variables with mean µ x 1 and µ x 2 , respectively. The number of equations required for the nth order polynomial is [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF]:

N eq = 2n i=1 (N r v -i + 1)! i!(N r v -1) (2.40)
where N eq is the number of the equations required for the nth order approximation with N r v random variables. The above equation shows that each sigma point adds N r v + 1 unknowns to each of the N eq equations. Consequently, the number of sigma points m is the next integer to the ratio between N eq and N r v + 1:

m = N eq N r v + 1 (2.41)
From the above equations, the number of equations N eq increases rapidly with N r v .

Higher-order UT

It is obvious that the accuracy of the UT method can be improved by taking more terms into account from the Taylor series approximation in (2.7). For example, a fourth-order approximation in the case of one random variable in UT calculation (one-dimensional UT) is [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF]:

h(x) = a 1 x + a 2 x2 + a 3 x3 + a 4 x4 (2.42)
Equation (2.39) shows that a fourth-order approximation requires eight equations to find the weights and sigma points. If the random variable is a normally distributed random variable, the eight equations are:

E[x] = m-1 i=1 w i S i = 0 E[x 2 ] = m-1 i=1 w i S 2 i = σ 2 x E[x 3 ] = m-1 i=1 w i S 3 i = 0 E[x 4 ] = m-1 i=1 w i S 4 i = 3σ 4 x E[x 5 ] = m-1 i=1 w i S 5 i = 0 E[x 6 ] = m-1 i=1 w i S 6 i = 15σ 6 x E[x 7 ] = m-1 i=1 w i S 7 i = 0 E[x 8 ] = m-1 i=1 w i S 8 i = 105σ 8 x (2.43)

Uncertainty quantification using the Unscented Transformation (UT)

There are four weights and four sigma points required to solve the eight equations which cannot be solved directly. It requires using an alternative method to simplify the calculation of the fourth-order UT's sigma points and weights. In [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF], the author proposes using the relationship between the UT technique and the integral calculation of moment to solve the equations. The idea is using the Gaussian Quadrature method [12] to approximate the integral by selecting the optimal abscissas at which to evaluate the function, and it is given by:

+∞ -∞ f (x)W (x)d x ≈ m-1 i=1 w i f (S i ) = E[ f (x)] (2.44)
where f (x) is the function, W (x) is the weighting function representing the continuous PDF, w i are the weights and S i are the optimal abscissas which correspond with the sigma points in UT. Orthogonal polynomials ψ n (x) are a class of polynomials defined in a closed interval [x 1 , x 2 ] and are related to the weighting function W (x) by the equation

x 2 x 1 W (x)ψ i (x)ψ j (x)d x = 0 i ̸ = j, i, j = 0, 1, 2 . . . (2.45)
Then the zeros of ψ n (x) of degree n are the optimal abscissas and correspond to the sigma points S i in UT, and the number of sigma points S i is equal to n. ψ n (x) is defined as follows:

ψ n (x) = 1 c n W (x) d n d x n [W (x)G(x) n ] (2.46)
where c n is a constant and G(x) is a polynomial independent of n. The weights w i are calculated as follows:

w i = 1 dψ n d x | x=S i x 2 x 1 ψ n (x) x -S i W (x)d x (2.47)
Choosing which kind of polynomial ψ n (x) depends on the W (x). For example, if W (x) is a Gaussian distribution, then ψ n (x) is a Hermite polynomial. If W (x) is an exponential distribution, a Laguerre polynomial is chosen. The form of the Hermite polynomial is:

ψ n (x) = (-1) n e x 2 2 d n d x n [e -x 2 2 ]
(2.48)

ψ n+1 (x) = xψ n (x) -nψ n-1 (x) (2.49)
The sigma points and weights for the fourth-order UT obtained using Hermite polynomials at n = 5 within the range (-∞, +∞) are given by:

S 0 = 0 S 1 = 5 -10σ x S 2 = -5 -10σ x S 3 = 5 + 10σ x S 4 = -5 + 10σ x w 0 = 1 -w 1 -w 2 -w 3 -w 4 w 1 = w 2 = 3 4(5 -10) 2 w 3 = w 4 = 3 4(5 + 10) 2 (2.50)
If the PDF of the random variable is an exponential distribution, Laguerre polynomials are used to find the parameters in UT. The generating function for a Laguerre polynomial is given by [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]:

L n (x) = e n n! d n d x n [x n e -x ] (2.51)
The recurrence relation for L n (x) is:

(n + 1)L n+1 (x) = (2n + 1 -x)L n (x) -nL n-1 (x) (2.52)
The sigma points and weights of fourth-order UT for an exponential distribution are found from Laguerre polynomial at n = 5 within the range (0, +∞) [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF]:

S 0 = 0.596425771σ x S 1 = 0.2635603197σ x S 2 = 1.413403059σ x S 3 = 2.64080084σ x S 4 = 7.085810006σ x w 0 = 0.7594244976 × 10 -1
w 1 = 0.5217556100 w 2 = 0.3986668121

w 3 = 0.3611758646 × 10 -2 w 4 = 0.2336997028 × 10 -4 (2.53)
Uncertainty quantification using the Unscented Transformation (UT)

The higher the order of an approximation, the more the weights and sigma points needed to be calculated satisfying the moment equations [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF].

Results and discussions

Here in the Microwave Studio based on FIT from CST is used for the numerical simulations. A commercially available dual-band mobile phone -900 MHz and 1800 MHz frequency bands-is investigated, as presented in Fig.

Example 1: using second-order UT for UQ with one random variable in the calculation

The results obtained by second-order UT correspond to three runs for the case of only one random variable in the calculation. In the case of random variable x 1 in the calculation, the results for 900MHz and 1800MHz are presented in Fig. 2.3 and Fig. 2.4, respectively. At 900MHz, the result of UT is 0.895 ± 0.73%, and the result of MCM corresponding to 1000 runs is 0.895 ± 0.74%. At 1800MHz, the result of UT is 1.719 ± 1.51%, and the results of MCM corresponding to 1000 runs is 1.718 ± 1.47%. When the random variable in the calculations is x 2 , the results obtained by MCM and second-order UT are given in Table 2.1. The results of Fig. 2.3, Fig. 2.4 and Table 2.1 show that MCM and second-order UT give similar results in the case of only one random variable in the calculation.

Example 2: using second-order UT for UQ with two random variables in the calculation

In the case of two random variables, x 1 and x 2 in the calculation, the results of 900MHz and 1800MHz obtained by MCM and UT are presented in Fig. 2.5 and Fig. 2.6, respectively. At 900MHz, the result of UT corresponding to 36 runs is 0.896 ± 1.29%, and the result of MCM corresponding to 1000 runs is 0.895 ± 1.26%. At 1800MHz, the result of UT in the case of two random variables in the calculation is 1.719 = ±1.71%, and the result of MCM running 1000 times is 1.719 ± 1.98%. The results whether for 900MHz or 1800MHz obtained by MCM and second-order UT are similar.

Example 3: using fourth-order UT for UQ with one random variable in the calculation

For fourth-order UT, the sigma points are calculated by the method proposed in [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF], and the sigma points are given by equation (2.44) using this method. When there is only one random variable x 1 in the calculation, the mean and standard variance of the output for 900MHz and 1800MHz calculated by MCM, second-order UT and fourthorder UT are shown in Fig. 2.7 and Fig. 2.8, respectively. The result of fourth-order UT is 0.895 ± 1.35% at 900MHz and 1.719 ± 3.57% at 1800MHz. The results displayed show that the mean of the output calculated by second-order UT and fourth-order UT are similar for both 900MHz and 1800MHz, however, the standard variance of the output calculated by second-order UT is closer to the results of MCM than the standard deviation calculated by fourth-order UT. Theoretically, the fourth-order UT should be more accurate than the second-order UT since higher-order of Taylor series are considered. However, since the equations of fourth-order UT cannot be calculated directly, the sigma points of fourth-order UT are approximated by the method proposed in [START_REF] Ajayi | Direct computation of statistical variations in electromagnetic problems[END_REF], and some error is existing. The error results in the poor performance of fourthorder UT. In contrast, the sigma points of the second-order UT can be obtained directly by solving the equations without approximation or error. Fig. 2.8 In the case of one random variable, the mean and standard deviation of the output calculated by MCM, second-order UT and fourth-order UT, respectively, for 1800MHz

Chapter 3

Uncertainty quantification using non-intrusive methods

Introduction

In this chapter, two other non-intrusive UQ methods are introduced. These two methods are the SC method and NIPC. SC expansion and PC expansion are very similar methods because both of them have the ability to produce functional representations of stochastic variability. SC forms interpolation functions for known coefficients, and requires the use of structured collocation point sets derived from tensor product or sparse grids [START_REF] Eldred | Evaluation of non-intrusive approaches for Wiener-Askey generalized polynomial chaos[END_REF]. PC expansion, on the other hand, estimates coefficients for known orthogonal polynomial basis functions based on a set of response function evaluations, using sampling, linear regression, tensor product quadrature, or Smolyak sparse grid approaches. The performances of UT, SC and NIPC for UQ in SAR calculation using a commercial available mobile phone model are compared and analysed. To quantify the relative importance of the uncertain inputs in SAR calculation, one of variance-based SA methods related to PC expansion is used.

The formulations of multidimensional UQ methods are complex, and the number of simulations required for UQ increases with the number of input parameters. However, in SAR calculation there are usually several random variables due to the lack of precise knowledge of the material or geometry of the components of the mobile phone being modelled. To solve this problem, combined uncertainties is investigated to be used in SAR calculation. Combining uncertainties is an approach following the idea of the MoM given in chapter 1.

Stochastic collocation (SC) method

Lagrange polynomials [START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF]35] interpolate a set of points in a single dimension using the functional form :

L i (x) = k j̸ =i x -x j x i -x j (3.1)
where it is evident that L i is 1 at x = x i , and is 0 for each of the points x = x j .

SC expansion [START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF] is formed as a sum of a set of multidimensional Lagrange interpolation polynomials, and the function f (x) can then be expressed as follows:

f (x) = k i=1 f (x i )L i (x) (3.2)
where L i are the Lagrange interpolation polynomials, x i are so-called particular collocation points which need to be chosen based on a certain rule, and k is the number of collocation points. Since these polynomials have the feature of being equal to 1 at their particular collocation point and 0 at all other points, the coefficients of the expansion are simply the response values at each of the collocation points. For interpolation in multiple dimensions, a tensor-product approach [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] is used and the form is:

f (x) = m i 1 j 1 =1 • • • m i N j N =1 f (x i 1 j 1 , . . . , x i N j N )(L i 1 j 1 ⊗ . . . ⊗L i N j N ) (3.3)
where multi-index {i 1 , i 2 , ..., i N } T ∈Z + N and x = {x i , . . . , x N } For interpolation in one dimension N = 1, the mean µ y of the output y is then formed as follows:

µ y = k i=1 f (x i ) L i (x)ρ(x)d x (3.4)
where ρ(x) is the probability density function of the random variable x. These polynomials L i have the feature of being equal to 1 at their particular collocation point and 0 at all other points. µ y can be estimated by numerical integration or quadrature as follows:

µ y = 〈 f (x)〉 ≈ k i=1 f (x i )w i (3.5)
where 〈 〉 denotes the inner product, w i and x i are the quadrature weights and points, respectively. Similarly, the variance σ 2 y of y can be derived for stochastic collocation:

σ 2 y = 〈 f 2 (x)〉 -µ 2 y ≈ k i=1 k j=1 f (x i ) f (x j )〈L i (x)L j (x)〉 -µ 2 y = k i=1 f 2 (x i )w i -µ 2 y (3.6)
The collocation points x i are chosen by quadrature rule. When the dimension of the random variables N > 1 (i.e., multiple random variables), the mean µ y and variance σ 2 y of the multidimensional SC expansion are estimated by a full tensor product quadrature.

The formula of the full tensor product quadrature are as follows:

µ y = m i 1 j 1 =1
...

m i N j N =1 f (x i 1 j 1 , ..., x i N j N )(w i 1 j 1 ⊗ ...⊗w i N j N ) (3.7) σ 2 y = m i 1 j 1 =1
...

m i N j N =1 f 2 (x i 1 j 1 , ..., x i N j N )(w i 1 j 1 ⊗ ...⊗w i N j N ) -µ 2 y (3.8)
where multi-index {i 1 , i 2 , ..., i N } T ∈Z + N . It is obvious that the above product needs N n=1 m i n function evaluations. The number of collocation points in the tensor grid grows exponentially with the number of input random variables. Using full tensor product quadrature method for faces the curse of dimensionality.

Polynomial chaos (PC) expansion

The PhD thesis of Edwards [21] gives detailed information about PC expansion. The concept of Homogeneous Chaos was first introduced by Wiener in 1938 [START_REF] Wiener | The Homogeneous Chaos[END_REF]. Ghanom and Spanos [START_REF] Ghanem | Stochastic Finite Elements with multiple random non-Gaussian Properties[END_REF][START_REF] Ghanem | Ingredients for a General Purpose Stochastic Finite Element Formulation[END_REF][START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] pioneered PC expansion and applied it to problems in mechanics. In the generalised PC expansion the output of interest is approximated as a function of certain orthogonal basis polynomials. The choice of the basis polynomials is dependent on the probability distributions of the uncertain input parameters [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF]. It is possible to solve stochastic differential equations, involving the approximate function, by using the orthogonality of the polynomial basis set. The orthogonality of the basis polynomials reduces the stochastic differential equations to a set of deterministic differential equations, which can be solved numerically. Xiu and Karniadakis [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF] found this approach to be computationally cheaper than the Monte Carlo Method. They found the mean and variance for a specific problem 500 times faster than the Monte Carlo Method, for Uncertainty quantification using non-intrusive methods the same accuracy. They noted however that the methods' efficiency is problem specific [21].

Homogeneous chaos

Wiener [START_REF] Wiener | The Homogeneous Chaos[END_REF] proposed the use of a homogeneous chaos expansion to represent a process, which depends on normally distributed uncertain input parameters, in terms of Hermite polynomials. Cameron and Martin [START_REF] Cameron | The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals[END_REF] proved that the expansion can be used to approximate any second-order linear function. This means that the homogeneous expansion can be used to expand second-order random processes in terms of Hermite polynomials. Second-order random processes are processes with finite variance, which encompasses most physical processes. A general second-order random process M (θ ) which depends on a normally distributed random event θ can be represented as the homogeneous chaotic expansion of M (θ ) as

M (θ ) = a 0 H 0 + ∞ i 1 =1 a i 1 H 1 (x i 1 (θ )) + ∞ i 1 =1 i 1 i 2 =1 a i 1 i 2 H 2 (x i 1 (θ ), x i 2 (θ )) + ∞ i 1 =1 i 1 i 2 =1 a i 1 i 2 i 2 i 3 =1 a i 1 i 2 i 3 H 3 (x i 1 (θ ), x i 2 (θ ), x i 3 (θ )) + . . . (3.9)
where H n (x i 1 (θ ), . . . , x i n (θ )) is the Hermite polynomial chaos of order n in terms of the n standard normally distributed variables (x i 1 (θ ), . . . , x i n (θ )) which have zero mean and unit variance. When x = (x i 1 (θ ), . . . , x i n (θ )), the Hermite polynomials are defined as

H n (x) = e 1 2 x T x (-1) n ∂ n ∂ x i 1 . . . ∂ x i n e -1 2 x T x (3.10)
For convenient notation, equation (3.9) can be rewritten as

M (θ ) = ∞ i=0 a i ψ i (x) (3.11)
where there is a one-to-one correspondence between H i and ψ i , and a i are constant coefficients of the expansion to be determined.

The polynomial basis ψ i (x) is a set of orthogonal polynomials which satisfy the orthogonality condition 〈ψ i (x)ψ j (x)〉 = 〈ψ 2 i 〉δ i j (3.12) where δ i j is the Kronecker delta [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF], and 〈〉 denotes the inner product. The inner product operation defined in the Hilbert space of the random vector x

〈 f (x)g(x)〉 = f (x)g(x)W (x)d x (3.13)
where W (x) is the PDF of the standard normally distributed variables corresponding to the Hermite polynomial basis set. It is given by:

W (x) = 1 (2π) n 1 2 x T x (3.14)
Homogeneous chaos was first applied by Ghanem and Spanos [START_REF] Ghanem | Stochastic Finite Elements with multiple random non-Gaussian Properties[END_REF][START_REF] Ghanem | Ingredients for a General Purpose Stochastic Finite Element Formulation[END_REF][START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. They applied this theory to problems in mechanics when the random variables are Normally distributed random variables. Xiu and Karniadakis [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF] extended this approach to a more general case so as the expansion can accurately represent stochastic processes that depend on random variables, which are not necessarily Normally distributed. It is called the Wiener-Askey Chaos expansion or the generalized Polynomial Chaos (gPC) expansion, which is based on polynomial bases from the Askey scheme.

The generalized polynomial chaos (gPC) method

The gPC method proposed by Xiu and Karniadakis can employ more types of orthogonal polynomials based on the distribution of random variables. It means that any secondorder random process M (θ ) depending on the random event θ can be represented as

M (θ ) = ∞ i=0 a i ψ i (x(θ )) (3.15)
As indicated already, this is the Wiener-Askey Chaos expansion when the polynomial basis set ψ i corresponding to the random variables x = (x i 1 (θ ), . . . , x i n (θ )) is chosen from the Askey-scheme as shown in Table 3.1. The polynomial basis sets ψ i are all orthogonal, which satisfies equation (3.12).

Using non-intrusive polynomial chaos for UQ

An approach called NIPC similar to SC is proposed in [START_REF] Eldred | Evaluation of non-intrusive approaches for Wiener-Askey generalized polynomial chaos[END_REF]. NIPC methods allow the use of simulations as black boxes for UQ studies, and involve the calculation of chaos expansion coefficients based on a set of response function evaluations. Based on the idea that a set of polynomials is used as an orthogonal basis to approximate the functional form between the stochastic response output and each of its random inputs, gPC expansion can be simplified as

y = f (x) ≈ P i=0 a i ψ i (x) (3.16)
where f is a mathematical model of a physical system depending on input parameters x = {x 1 , x 2 , ..., x N } T ∈ R N , y is the output of the model, and a i are polynomial coefficients. Here, the infinite sum has been truncated at P, and P + 1 is the total number of terms given by:

P + 1 = (r + N )! r!N ! (3.17)
where r is the order of PC and N is the number of the random variables. For example, when the random vector x contains normal distribution random variables, ψ i are Hermite polynomials. When the dimension of x is N = 1, the random variable is x = {x 1 }, and the order of PC expansion is r = 2, then the total number of terms is P + 1 = 3, and ψ i (x) are:

ψ 0 (x 1 ) = 1 ψ 1 (x 1 ) = x 1 ψ 2 (x 1 ) = x 1 2 -1 (3.18)

Polynomial chaos (PC) expansion

47

When there are multiple independent random variables i.e., x = {x 1 , x 2 , . . . , x N } T ∈ R N , the multidimensional gPC basis ψ i is constructed by tensor products of the corresponding one-dimensional polynomial:

ψ s (x 1 , x 2 , ..., x N ) = N i=1 ψ s i (x i ) s = {s i , i = 1, ..., N }, s i ≥ 0, N i=1 s i ≤ r (3.19)
where ψ s i (x i ) denote one-dimensional polynomials. For example, N = 2, r = 2, and the probability distributions of the random variables are normal distributions, the two-dimensional Hermite polynomials ψ s (x 1 , x 2 ) are:

ψ 00 (x 1 , x 2 ) = ψ 0 (x 1 )ψ 0 (x 2 ) = 1 ψ 10 (x 1 , x 2 ) = ψ 1 (x 1 )ψ 0 (x 2 ) = x 1 ψ 01 (x 1 , x 2 ) = ψ 0 (x 1 )ψ 1 (x 2 ) = x 2 ψ 20 (x 1 , x 2 ) = ψ 2 (x 1 )ψ 0 (x 2 ) = x 2 1 -1 ψ 11 (x 1 , x 2 ) = ψ 1 (x 1 )ψ 1 (x 2 ) = x 1 x 2 ψ 02 (x 1 , x 2 ) = ψ 0 (x 1 )ψ 2 (x 2 ) = x 2 2 -1 (3.20)
where the total number of terms is P + 1 = 6. Another example, N = 3, r = 2 and the three-dimensional Hermite polynomials ψ s (x 1 , x 2 , x 3 ) are:

ψ 000 (x 1 , x 2 , x 3 ) = ψ 0 (x 1 )ψ 0 (x 2 )ψ x 3 (x 3 ) = 1 ψ 100 (x 1 , x 2 ) = ψ 1 (x 1 )ψ 0 (x 2 )ψ 0 (x 3 ) = x 1 ψ 010 (x 1 , x 2 , x 3 ) = ψ 0 (x 1 )ψ 1 (x 2 )ψ 0 (x 3 ) = x 2 ψ 001 (x 1 , x 2 , x 3 ) = ψ 0 (x 1 )ψ 0 (x 2 )ψ 1 (x 3 ) = x 3 ψ 110 (x 1 , x 2 , x 3 ) = ψ 1 (x 1 )ψ 1 (x 2 )ψ 0 (x 3 ) = x 1 x 2 ψ 101 (x 1 , x 2 , x 3 ) = ψ 1 (x 1 )ψ 0 (x 2 )ψ 1 (x 3 ) = x 1 x 3 ψ 011 (x 1 , x 2 , x 3 ) = ψ 0 (x 1 )ψ 1 (x 2 )ψ 1 (x 3 ) = x 2 x 3 ψ 200 (x 1 , x 2 , x 3 ) = ψ 2 (x 1 )ψ 0 (x 2 )ψ 0 (x 3 ) = x 2 1 -1 ψ 020 (x 1 , x 2 , x 3 ) = ψ 0 (x 1 )ψ 2 (x 2 )ψ 0 (x 3 ) = x 2 2 -1 ψ 002 (x 1 , x 2 , x 3 ) = ψ 0 (x 1 )ψ 0 (x 2 )ψ 2 (x 3 ) = x 2 3 -1 (3.21)
where the total number of terms is P+1 = 10. Using the orthogonality of the polynomials, the statistics of the output y such as mean µ y and variance σ 2 y can be calculated as:

µ y = 〈 f (x), 1〉 = P i=0 a i 〈ψ i , 1〉 = P i=0 a i 〈ψ i , ψ 0 〉 = P i=0 a i δ i0 = a 0 (3.22)
The variance σ 2 y is calculated as [START_REF] Eldred | Evaluation of non-intrusive approaches for Wiener-Askey generalized polynomial chaos[END_REF]:

σ 2 y = 〈 f (x), f (x)〉 -〈 f (x), 1〉 2 = P i=0 P j=0 a i a j 〈ψ i , ψ j 〉 -a 0 2 = P i=0 P j=0 a i a j δ i j 〈ψ 2 i 〉 -a 0 2 = P i=1 a 2 i 〈ψ 2 i 〉 (3.23)
When the coefficients of polynomials are calculated, then the statistics of output are formed by the above equations. Therefore, the important step for NIPC is to calculate each coefficient a i .

PC is expanded using orthogonal polynomials which are functions of independent standard random variables x. If the original random variables ξ are not independent standard random variables, they should be transformed to the proper form. In [START_REF] Eldred | Evaluation of non-intrusive approaches for Wiener-Askey generalized polynomial chaos[END_REF], methods for transforming the original random variables ξ to independent standard random variables x are presented. For problems directly involving independent normal, uniform, exponential, beta, and gamma distributions for input random variables, conversion to standard form involves a simple linear scaling transformation (to the form of the density functions in Table 3.1) and then the corresponding chaos can be employed. For correlated normal, uniform, exponential, beta, and gamma distributions, the same linear scaling transformation is applied, followed by application of the inverse Cholesky factor of the correlation matrix.

For other distributions with a close explicit relationship to variables supported in the Askey scheme, a non-linear transformation is employed to transform to the corresponding Askey distributions, i.e., normals and uniforms and the corresponding chaos polynomials, i.e., Hermite and Legendre are employed. For other less directly-related distributions, the non-linear Nataf transformation is employed to transform to uncorrelated standard normals as described below, and Hermite chaos is employed. The transformation from correlated non-normal distributions to uncorrelated standard normal distributions is denoted as x = T (ξ), with the reverse transformation denoted as ξ = T -1 (x). These transformations are non-linear in general, and possible approaches include the Rosenblatt [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF], Nataf [START_REF] Der Kiureghian | Structural reliability under incomplete probability information[END_REF], and Box-Cox transformations [START_REF] Box | An analysis of transformations[END_REF]. For simplicity, it is assumed from now on that the components of the original input random variable ξ have been already scaled, i.e., x = ξ. All input parameters are assumed independent standard random variables.

To obtain the statistics of the output of the model by NIPC, as mentioned above, the key issue is to calculate the constant coefficients a i for which there are two main methods: Galerkin projection and the regression method. The Galerkin projection or projection method has been introduced in Ghanem and Ghiocel (1998) [START_REF] Ghanem | Stochastic seismic soil-structure interaction using the homogeneous chaos expansion[END_REF], [START_REF] Ghiocel | Stochastic finite element analysis of seismic soilstructure interaction[END_REF] [START_REF] Ghiocel | Stochastic finite element analysis of seismic soilstructure interaction[END_REF] in the context of seismic soil-structure interaction, Ghanem et al. (2000) [START_REF] Ghanem | Modal properties of a space-frame with localized system uncertainties[END_REF] for the study of the eigenmodes of a spatial frame structure and Fieldet al. (2000) [START_REF] Field | A non deterministic shock and vibration application using polynomial chaos expansion[END_REF], [START_REF] Field | Numerical methods to estimate the coefficients of the polynomial chaos expansion[END_REF] [START_REF] Field | Numerical methods to estimate the coefficients of the polynomial chaos expansion[END_REF] in a non linear shock and vibration analysis. Le Maître et al. (2002) [61] make use of the same so-called non intrusive spectral projection (NISP) to solve a natural convection problem in a square cavity [12]. The Galerkin projection is given as follows:

〈 f (x), ψ k 〉 = 〈 P i=0 a i ψ i , ψ k 〉 (3.24)
employing the orthogonality of polynomial basis ψ i :

〈 f (x), ψ k 〉 = a k 〈ψ 2 k 〉 (3.25) a k = 〈 f (x), ψ k 〉 〈ψ 2 k 〉 (3.26)
the inner product can be estimated by quadrature method or MCM. For example, if the distributions of the random variables are Gaussian, then Gauss-Hermite quadrature points can be used to estimate the inner products. When the dimension of random variable N > 1, the full tensor product quadrature method as we mentioned in the SC method should be applied here also. When there are multiple uncertain inputs, using the full tensor product quadrature method to calculate faces the problem known as the curse of dimensionality [START_REF] Walters | Uncertainty analysis for fluid dynamics with applications[END_REF]. An alternative, the so-called sparse quadrature scheme originally introduced by Smolyak (1963) [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF], has been proposed to overcome this curse.

The regression method is based on the least square minimization of the error between the model output and its approximation. [START_REF] Tatang | An efficient method for parametric uncertainty analysis of numerical geophysical models[END_REF] and Isukapalli (1999) [START_REF] Isukapalli | Uncertainty Analysis of Transport-Transformation Models[END_REF] under the name "probabilistic collocation methods" [12]. The least squares problem is: given a matrix Ψ ∈ R N ×(P+1) , N > (P + 1) and a vector y ∈ R N , find a vector a ∈ R (P+1) such that Ψa is the best approximation to y in the least squares sense. Mathematically speaking, we wish to solve the following minimization problem:

min = ∥ Ψa -y ∥ 2 (3.27)
It is assumed that the rank of matrix Ψ is equal to P + 1. Several algorithms can be applied to compute the least squares solution.

Equation (3.16) can be rewritten as follows:

y = Ψa (3.28)
where y = ( f (x 1 ), . . . , f (x N )) T , Ψ = (ψ j (x i ), i = 1, . . . , N , j = 0, . . . , P) ∈ R N ×(P+1) is the matrix of polynomials, N is the number of sampling points, and the vector of the coefficients is a = (a 0 , . . . , a P ) T . It is calculated as a least-square problem with solution as follows:

a = (Ψ T Ψ) -1 Ψ T y (3.29)
The above equation is only valid for a full-rank information matrix. A necessary condition is that the size N of the experimental design is not less than the number P + 1 of PC coefficients to estimate. The disadvantage is that it requires a dramatically increasing number of calculations when the number of uncertain input parameters in the model increases. To overcome this difficulty, sparse least squares schemes may be employed [5]. These methods are aimed at identifying a small set of significant coefficients in the model response approximation. Eventually, a sparse polynomial response surface, i.e. a

Polynomial chaos (PC) expansion

51

polynomial which only contains a small number of non-zero coefficients, is obtained by means of a reduced number of possibly costly model evaluations.

Sensitivity analysis (SA) using PC expansion

The modelling of complex systems usually requires a large number of input parameters such as the permittivity and conductivity of the material of a mobile phone. Defining a proper probabilistic model for these parameters is a complex work. However, in most real-world problems, only a limited number of input parameters happens to influence the response randomness significantly [12]. SA of a model aims at quantifying the relative importance of each input parameter. Methods of SA are usually classified into two categories [12]:

• Local sensitivity analysis concentrates on the local impact of input parameters on the model. It is based on the computation of the gradient of the response with respect to its parameters around a nominal value. Numerous techniques have been developed to compute the gradient efficiently, including finite-difference schemes, direct differentiation or adjoin differentiation methods [START_REF] Cacuci | Sensitivity and Uncertainty Analysis: Theory[END_REF].

• Global SA tries to quantify the output uncertainty due to the uncertainty in the input parameters, which are taken singly or in combination with others.

Many papers have been devoted to global SA in the last two decades or so. A good state-of-the-art of the techniques is available in [START_REF] Saltelli | Sensitivity Analysis[END_REF] [69], which gathers the methods into two groups [12]:

• Regression-based methods: the Standardized Regression Coefficients (SRC) are based on a linear regression of the output on the input vector. The input/output Pearson correlation coefficients measure the effect of each input variable by the correlation it has with the model output. The Partial Correlation Coefficients (PCC) are based on results of regressions of the model on all input variables except one. These coefficients are useful to measure the effect of the input variables if the model is linear, i.e., if the coefficient of determination R 2 of the regression is close to one. In case of non-linearity, they fail to properly represent the response sensitivities. However, in case of monotonicity of the model with respect to the input parameters, the rank transform can be used, leading to the so-called standardized rank regression and partial rank correlation coefficients. As a whole, in case of general non-linear non-monotonic models, these approaches fail to produce satisfactory sensitivity measures (Saltelli and Sobol', 1995) [START_REF] Saltelli | About the use of rank transformation in sensitivity of model output[END_REF].

• Variance-based methods: these methods aim at decomposing the variance of the output as a sum of contributions of each input variable, or combinations thereof. They are sometimes called ANOVA techniques for "Analysis Of Variance". The correlation ratios proposed in McKay (1995) [START_REF] Mckay | Evaluating prediction uncertainty[END_REF] enter into this category. They are formulated as conditional variances and usually evaluated by crude Monte Carlo simulation or Latin Hypercube Sampling. The Fourier amplitude sensitivity test (FAST) indices (Cukier et al., 1978 [72]; Saltelli et al., 1999 [73]) and the Sobol' indices (Sobol', 1993 [74]; Saltelli and Sobol', 1995 [START_REF] Saltelli | About the use of rank transformation in sensitivity of model output[END_REF]; Archer et al., 1997 [76]), see also the review in Sobol' and Kucherenko (2005) [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models[END_REF], are intended to represent the sensitivities for general models. The Sobol' indices are in practice computed using the MCM, which makes them hardly applicable for computationally demanding models, e.g., finite element models in engineering mechanics.

Sobol decomposition [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] is a commonly used variance-based SA method. Its sensitivity indices, i.e., Sobol' indices are generally calculated by MCM. In [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF], Sudret proposes to use PC expansion to calculate the Sobol' indices for SA. The so-called Sobol' decomposition of a function f (x 1 , . . . , x N ) depending on inputs x = {x 1 , x 2 , . . . , x N } T ∈ R N is:

f (x 1 , . . . , x N ) = f 0 + N i=1 f i (x i ) + 1≤i< j≤N f i j (x i , x j ) + • • • + f 1,2,...,N (x 1 , . . . , x N ) (3.30)
where f 0 is a constant that corresponds to the mean of the output and each summand f i 1 ,...,i m (x i 1 , . . . , x i m ) over any of its independent variables is zero, i.e.:

1 0 f i 1 ,...,i m (x i 1 , . . . , x i m )d x i k = 0, 1 ≤ k ≤ m. (3.31)
The properties of Sobol decomposition are as follows [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]:

• The sum in equation (3.30) contains a number of summands equal to

N i=1 N i = N i=1 (N )! (N -i)!i! = 2 N -1 (3.32)

Polynomial chaos (PC) expansion

53

• The constant f 0 is the mean value of the function:

f 0 = f (x)dx (3.33)
where dx stands for d x 1 . . . d x N for the sake of simplicity.

• The summands are orthogonal to each other as follows:

f i 1 ,...,i m (x i 1 , . . . , x i m ) f k 1 ,...,k t (x k 1 , . . . , x k t )dx = 0, (i 1 , . . . , i m ) ̸ = (k 1 , . . . , k t ) (3.34)
Considering that the input parameters are independent random variables, the total variance D of the model response y = f (x) is given by: where the partial variances are :

D = σ 2 y = f 2 (x)d x -f 2 0 ( 3
D i 1 ,...,i m = f 2 i 1 ,...,i m (x i 1 , . . . , x i m )d x i 1 . . . d x i m 1 ≤ i 1 • • • < i m ≤ N , m = 1, . . . , N (3.37)
The Sobol' indices are defined as follows:

SU i 1 ,...,i m = D i 1 ,...,i m /D (3.38)
Considering equation (3.36), these indices satisfy:

N i=1 SU i + 1≤i< j≤N SU i j + ...SU 1,2,...,N = 1 (3.39)
Each index SU i 1 ,...,i m is a sensitivity measure describing which amount of the total variance is due to the uncertainties in the set of input parameters {x i 1 , . . . , x i m }. The first-order indices SU i provide the influence of each parameter taken alone whereas the higher-order indices account for possible mixed influence of various parameters [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF].

According to the idea mentioned above, the index SU 1,2,...,k can be evaluated by PC expansion as follows:

SU 1,2,...,k = s a 2 s E[ψ 2 s (x 1 , x 2 , ..., x k )] σ 2 y , 1 ≤ k ≤ N (3.40)
where the input random variables are x = {x 1 , x 2 , ...,

x N } T ∈ R N ,
s = {s i , i = 1, ..., N }, s i ≥ 0, N i=1 s i ≤ r
is larger than 0, so as ψ 00 or ψ 000 are not included. For example, N = 2 and r = 2, the sensitivity indices are

SU 1 = a 2 10 ψ 2 10 (x 1 , x 2 ) + a 2 20 ψ 2 20 (x 1 , x 2 ) σ 2 y (3.41) SU 2 = a 2 01 ψ 2 01 (x 1 , x 2 ) + a 2 02 ψ 2 02 (x 1 , x 2 ) σ 2 y (3.42) SU 12 = a 2 11 ψ 2 11 (x 1 , x 2 ) σ 2 y (3.43)
If N = 2 and r = 3, then the sensitivity indices are

SU 1 = a 2 10 ψ 2 10 (x 1 , x 2 ) + a 2 20 ψ 2 20 (x 1 , x 2 )) + a 2 30 ψ 2 30 (x 1 , x 2 ) σ 2 y (3.44) SU 2 = a 2 01 ψ 2 01 (x 1 , x 2 ) + a 2 02 ψ 2 02 (x 1 , x 2 ) + a 2 03 ψ 2 03 (x 1 , x 2 ) σ 2 y (3.45) SU 12 = a 2 11 ψ 2 11 (x 1 , x 2 ) + a 2 12 ψ 2 12 (x 1 , x 2 ) + a 2 21 ψ 2 21 (x 1 , x 2 ) σ 2 y (3.46)

Combined uncertainties

55

As third example, N = 3 and r = 2, the sensitivity indices are

SU 1 = a 2 100 ψ 2 100 (x 1 , x 2 , x 3 ) + a 2 200 ψ 2 200 (x 1 , x 2 , x 3 ) σ 2 y (3.47) SU 2 = a 2 010 ψ 2 010 (x 1 , x 2 , x 3 ) + a 2 020 ψ 2 020 (x 1 , x 2 , x 3 ) σ 2 y (3.48) SU 3 = a 2 001 ψ 2 001 (x 1 , x 2 , x 3 ) + a 2 002 ψ 2 002 (x 1 , x 2 , x 3 ) σ 2 y (3.49) SU 12 = a 2 110 ψ 2 110 (x 1 , x 2 , x 3 ) σ 2 y (3.50) SU 13 = a 2 101 ψ 2 101 (x 1 , x 2 , x 3 ) σ 2 y (3.51) SU 23 = a 2 011 ψ 2 011 (x 1 , x 2 , x 3 ) σ 2 y (3.52)
The Sobol' indices can be calculated by equation (3.40) when the polynomial coefficients are evaluated by the projection method or the regression method.

Combined uncertainties

The forms of multi-dimensional UQ methods are complicated, and the number of calculations required for UQ increases with the number of random variables. These problems exist in all UQ methods investigated in this thesis. Based on the guide to the expression of uncertainty in measurement, it may nevertheless be possible to combine the uncertainties induced by a given variable while the others are kept constant. This is possible when the variations in variables are small enough.

Following what has been said about the MoM in chapter 1, the steps of estimating the the first-order estimate of the variance are given next. The first-order estimate of the mean µ y1st and variance σ 2 y1st of the output of the model y = f (x) evaluated by the MoM are given by equation (1.28) and (1.31) as follows:

µ y1st = f (µ x ) σ 2 y1st = σ 2 x d f d x

56

Uncertainty quantification using non-intrusive methods where the dimension of input random vector is N = 1. When there are multiple input parameters x = (x 1 , . . . , x N ), N > 1, the first-order estimate of the mean is

µ y = f (µ 1 , . . . , µ N ) (3.53)
where µ y is mean of y, and µ i is the mean of x i .

From [START_REF]Joint Committee for Guides in Metrology, Evaluation of measurement data -Guide to the expression of uncertainty in measurement[END_REF], y = f (x 1 , x 2 , ..., x N ) with N input parameters x 1 , x 2 , ..., x N . Expansion of f around the expectation of the x i in a first-order Taylor series yields for small deviation of y around µ y in terms of small deviations of

x i around µ i y -µ y = N i=1 ∂ f ∂ x i (x i -µ i ) (3.54)
All higher-order terms are assumed to be negligible, and ( yµ y ) 2 is as follows:

( y -µ y ) 2 = N i=1 ( ∂ f ∂ x i ) 2 (x i -µ i ) 2 + (3.55) 2 N -1 i=1 N j=i+1 ∂ f ∂ x i ∂ f ∂ x j (x i -µ i )(x j -µ j )
The expectation of ( yµ y ) 2 is the variance σ 2 y of y and thus the σ 2 y of output y is as follows:

σ 2 y = N i=1 ( ∂ f ∂ x i ) 2 σ 2 i + 2 N -1 i=1 N j=i+1 ∂ f ∂ x i ∂ f ∂ x j σ i σ j ρ i j (3.56)
where σ i and σ j are uncertainties of x i and x j , respectively. ρ i j is the correlation coefficient between x i and x j . σ y is the first-order estimate of the uncertainty of y.

A small change △x i in input parameter x i induces the change in output y given by (△ y) i = (

∂ f ∂ x i )△x i . If the change in x i is △x i = σ i , (△ y) i = ( ∂ f ∂ x i )σ i . So σ 2
y can be estimated as follows:

σ 2 y ≈ N i=1 u 2 i ( y) + 2 N -1 i=1 N j=i+1 [u i ( y)u j ( y)ρ(x i , x j )] (3.57)
where u i ( y) and u j ( y) are the uncertainties of the output y induced by x i and x j , respectively. The correlation coefficient between x i and x j is -1 ≤ ρ(x i , x j ) ≤ +1.

If the estimates of them are independent, ρ(x i , x j ) = 0. If the ρ(x i , x j ) = +1, the maximum value of σ 2 can be obtained.

Results and discussions 57

Following a Taylor series expansion as mentioned above, the combined uncertainty u c or the overall uncertainty of the output may be written as follows:

u 2 c ≈ N i=1 u 2 i ( y) + 2 N -1 i=1 N j=i+1 [u i ( y)u j ( y)ρ(x i , x j )] (3.58)
where x = {x 1 , x 2 , ..., x N } T ∈ R N are the random variables, u i ( y) and u j ( y) are the uncertainties in output y induced by random variable x i and x j respectively, and ρ(x i , x j ) is the correlation coefficient between x i and x j . In the case that these variables are considered to be independent, the combined uncertainty u c is given by the first term of the summation, whereas in the general case the correlation between the variables should be accounted for.

In practice, even when the variables are apparently independent (e.g., the dielectric properties of the different materials), they may eventually be correlated, for example when the same instrument or measurement methodology has been applied to obtain the values. To evaluate the second term of the combined uncertainty, it becomes necessary to determine the correlation coefficient between each couple of random variables. Usually, little information is available to perform the evaluation. Nonetheless, an upper bound for the combined uncertainty may be derived when the correlation coefficient is assumed to be maximum in all possible configurations i.e., ρ(x i , x j ) = 1. The combined uncertainty equation is then expressed as follows:

u 2 c ≈ N i=1 u 2 i ( y) + 2 N -1 i=1 N j=i+1 [u i ( y)u j ( y)] (3.59)

Results and discussions

Uncertainty quantification using a simple mobile phone model

Before applying the non-intrusive methods introduced above for UQ in SAR calculation using a complex mobile phone model, first, the effectiveness of these UQ methods should be demonstrated. To demonstrate their effectiveness, the performances of these UQ methods should be compared with the traditional UQ method MCM. However, SAR calculation is time-consuming, and it is not practical to apply MCM for UQ in SAR calculation. A relatively simple example should be considered, thereby the MCM can be performed in the example.

For this purpose, these non-intrusive methods are used for UQ in the calculations which use the relatively simple mobile phone model presented in Fig. 2.1. The outputs are the two frequencies of the antenna for which the return losses are the minimum, as the examples presented in chapter 2. Since the mobile phone model is considered in free space configuration, and the model is relatively simple, the simulations are performed fast, and thus the MCM can be implemented. The results of UT, SC and NIPC are compared with the results obtained by MCM. The first random variable x 1 is the relative permittivity of the antenna support and the second random variable x 2 is the relative permittivity of the casing. Both of them are Gaussian probability distributions random variables as presented in chapter 2. The second-order UT is used, the number of collection points for SC is 3 in the case of one random variable in the calculation and 9 in the case of two uncertain input parameters in the calculation, and NIPC is divided into NIPC1 and NIPC2 depending on which method is used for calculating the polynomial coefficients. NIPC1 uses the regression method to calculate the coefficients, and NIPC2 uses the projection method to calculate the coefficients. The order of PC expansion is 2 for all cases. Fig. 3.1 gives the results obtained by MCM, second-order UT, SC and NIPC2 at 900MHz in the case of one random variable x 1 in the calculation. The result of MCM corresponding to 1000 runs is 0.895 ± 0.74%. The results of second-order UT, SC and NIPC2 are 0.895 ± 0.73%, 0.896 ± 0.89% and 0.895 ± 1.04%, respectively. At 1800MHz, the result obtained by MCM corresponding to 1000 runs is 1.718±1.47%. The results of second-order UT, SC and NIPC2 are 1.719 ± 1.51%, 1.720 ± 1.56% and 1.720 ± 2.14%, respectively, presented in Fig. 3.2. Fig. 3.3 and Fig. 3.4 show the results obtained by MCM and NIPC1 at 900MHz. There is one random variable x 1 in the calculations. The results of NIPC1 are obtained in two cases: without oversampling and with three times oversampling.

The numbers of sampling points or calculations for the two cases are 3 and 9, respectively. In the case of without oversampling, the results of NIPC1 at 900MHz using three different sets of sampling points are NIPC1 Case1=0.888 ± 0.80%, NIPC1 Case2=0.897 ± 2.51% and NIPC1 Case3=0.894 ± 2.06%, respectively. The results are given by Fig. 3.3. If three times oversampling is implemented in NIPC1, the results for different sampling sets are NIPC1 Case1=0.896 ± 0.23%, NIPC1 Case2=0.895 ± 0.88% and NIPC1 Case3=0.895 ± 0.46%, respectively, as presented in Fig. 3.4. It is obvious that implementing oversampling in NIPC1 improves its performance, making its results closer to the results of MCM (0.895 ± 0.74%). Without oversampling, the results of NIPC1 are not stable because of the influence of the bad sampling points. Oversampling

Results and discussions

59

is an approach to avoid or reduce the influence of bad sampling points. However, it increases the number of calculations and decreases the efficiency of NIPC1.

The results of two-dimensional UQ methods (with two uncertain input parameters in the calculation) at 900MHz and 1800MHz are given by Fig. 3.5 and Fig. 3.6, respectively. The results of MCM corresponding to 1000 runs is 0.895 ± 1.26%, and the results of second-order UT, SC and NIPC2 are 0.896 ± 1.29%, 0.896 ± 0.89%, 0.896 ± 1.79%, receptively, at 900MHz. Fig. 3.6 presents the results at 1800MHz in the same situation, the results of MCM corresponding to 1000 runs is 1.719 ± 1.98%, and the results of second-order UT, SC and NIPC2 are 1.719 ± 1.71%, 1.720 ± 2.17%, 1.721 ± 2.85%, respectively. Comparing with the results obtained by MCM, the results of second-order UT, SC, and NIPC2 for UQ are acceptable in the case of two random variables in the calculation.

With two uncertain parameters, the results obtained by two-dimensional NIPC1 and MCM are given in Fig. 3.7 and Fig. 3.8 for the two different frequencies, respectively. NIPC1 is performed in three different cases: without oversampling, with two times oversampling and with three times oversampling. At 900MHz, without oversampling, the result obtained by NIPC1 is 0.901 ± 0.42%, and the number of calculations is 6 for second-order PC expansion. With two times oversampling, the number of calculations required by NIPC1 is 12 and its result is 0.895 ± 0.94%. With three times oversampling, the number of calculations required by NIPC1 is 18, and the result is 0.895 ± 0.92%. At 1800 MHz, they are 1.706 ± 0.60%, 1.719 ± 1.76% and 1.719 ± 1.73%, respectively. With oversampling, the results of NIPC1 are approaching more to the results of MCM. These results once again prove that oversampling improves the performance of NIPC1, and it is necessary in NIPC1.

From the examples above, the conclusion is that comparing with the results obtained by MCM, the results for UQ using the simple mobile phone model obtained by secondorder UT, SC and NIPC2 are acceptable. The results obtained by UT are close to those of MCM comparing with SC and NIPC2. The performances of SC and NIPC2 are similar. The performance of NIPC1 associates with the number of sampling points, and oversampling improves its performance. Since all results are obtained in the case of variations in the input uncertain parameters not too large ( for x 1 is 16.67% and x 2 is 12.5%), the conclusion is drawn under this condition. In the following section, these non-intrusive UQ methods are applied in SAR calculation using a more complex mobile phone model.

Numerical model for SAR calculation

The model of a dual-band commercially available mobile phone called Neo FreeRunner [85] is used for SAR calculations in this section. Fig. 3.9 (a) shows the main components present in the CAD model of the dual-band commercially available mobile phone which was previously proposed as benchmark for the international inter-laboratory comparison. Table 3.2 provides the dielectric properties of the different components of the phone. Both numerical simulations and measurements of the mobile phone were undertaken by the laboratory leading in [10]. The critical components influencing the return loss of the antenna have been thoroughly examined, and some of the components found to have negligible electromagnetic contribution. For the purpose of the inter-laboratory comparison, some components present in the original CAD model have been discarded to avoid potential confusion in allotting material properties, and the so-called full CAD model consists of all the components given in Table 3.2 [10]. In this thesis the full CAD model is used for all SAR calculations. Fig. 3.9 (b) shows the mobile phone model placed against the right/cheek position of the SAM phantom. Table 3.

Uncertainty in SAR calculation

The SAR simulation process is modelled as y = f (x) with uncertain inputs x = {x 1 , x 2 , ..., x N } T ∈ R N and the outputs y, and the simulation process f is considered as a black box. The cases of one and two uncertain inputs in the calculations are considered, and the uncertain inputs are: (a) relative permittivity of the antenna support x 1 = 2.33 ± 10%, and (b) the relative permittivity of the bottom casing x 2 = 3 ± 10%. Both are assumed as Gaussian distributions and independent. When SAR calculations are implemented in Microwave Studio from CST, SAR can be normalized with respect to either accepted or stimulated power (SAR acc or SAR st i ). 3.10, which is in the planar cross section inside the phantom. Fig. 3.12 shows that for the variations in x 1 the SAR acc distributions are not much sensitive but the SAR acc distribution is relatively more sensitive at 1750 MHz than at 890 MHz.

The uncertainty of SAR acc can be quantified by the different UQ methods mentioned above. However, once the relationship between SAR acc and SAR st i is known, the uncer- tainty may be easily deduced from one another. The relationship between SAR acc and SAR st i is as follows:

SAR st i = (1 -S 2 11 )SAR acc (3.60)
where S 11 is the reflection coefficient.

Uncertainty quantification using non-intrusive methods

The uncertainties of SAR st i can be evaluated by combining uncertainties as mentioned in [START_REF]Joint Committee for Guides in Metrology, Evaluation of measurement data -Guide to the expression of uncertainty in measurement[END_REF] as follows:

u 2 SAR st i = ( ∂ SAR st i ∂ SAR acc ) 2 u 2 SAR acc + ( ∂ SAR st i ∂ S 11 ) 2 u 2 S 11 = (1 -S 2 11 ) 2 u SAR 2 acc + 4SAR 2 acc S 2 11 u 2 S 11 (3.61)
where u SAR st i , u SAR acc and u S 11 are the uncertainties of SAR st i , SAR acc and S 11 , respectively. The means and the uncertainties of S 11 and SAR acc need to be calculated by any of the UQ methods.

Here, an example about deducing the uncertainty of SAR st i from the uncertainty of SAR acc using SC method is given. SC can be replaced by any of the proposed UQ methods. Table 3.4 and Table 3.5 provide the SAR values from CST for these two cases respectively. SAR acc and SAR st i are normalized with respect to the same accepted powers and stimulated powers, respectively: 250 mW at 890 MHz and 125 mW at 1750 MHz. When the uncertain input x 1 changes, the reflection coefficients S 11 at 890 MHz are 0.4362, 0.5650 and 0.6651, and at 1750 MHz they are 0.5606, 0.2606 and 0.1976, respectively.

Based on the values provided in Table 3.4 and Table 3.5, u SAR st i are quantified by SC directly with 3 collocation points and combining uncertainties presented as equation [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF], respectively, and the results are given in Table 3.6 and Table 3.7. In order to be consistent, means and uncertainties of S 11 and SAR acc are evaluated by SC with 3 collocation points. At 890 MHz, the difference between the results of the two methods is acceptable. However, at 1750 MHz, u SAR st i quantified by the two approaches are different. In this case, the variation in S 11 is found to be important (about 55%) compared to only about 15% at 890 MHz. The combined uncertainties method is based on the assumption that the variations in the uncertain inputs are small and only in this case the high-order terms of Taylor series can be assumed negligible.

In the following sections, u SAR acc is quantified by the different UQ methods. However, u SAR st i can be evaluated by combining uncertainties under certain condition as mentioned in the example without implementing the same calculation twice. Table 3.8 presents the results of the mean and uncertainty of output-peak 10g average. The uncertainties in output induced by x 1 and x 2 taken independently are calculated. Second-order UT is used and 3 simulations are required. For SC, Gaussian quadrature is used to choose the 3 collocation points resulting in 3 simulations. The order of PC is r = 2. When the projection method is implemented to obtain PC coefficients, Gaussian quadrature is used with 3 points. The data in Table 3.8 shows that the mean and uncertainty of the output quantified by UT, SC and NIPC2 are in good agreement with only one uncertain input in the calculation. When linear regression method used for estimating PC coefficients represented as NIPC1, the order of PC is r = 2 and the sampling method used is LHS. Table 3.9 gives the results of NIPC1, and the random variable is x 1 . Oversampling is implemented, and the number of sampling points is n = 2(P + 1) , n = 3(P + 1), n = 4(P + 1) , n = 5(P + 1) and n = 6(P + 1), respectively. The results of Table 3.9 show that when increasing the number of sampling points, the results of NIPC1 are closer to the results of the other three methods. The performance of using the ordinary least-square estimator to evaluate PC coefficients is not good. Increasing the number of random samples seems to improve its performance. However, its efficiency is decreased.

Table 3.10 presents the results of calculation in the case of two random variables in SAR calculation. The two random variables x 1 and x 2 are assumed to be independent. The method presented in [START_REF] De Menezes | Efficient computation of stochastic electromagnetic problems using unscented transforms[END_REF] is used to calculate sigma points in UT and 36 simulations are required. For SC, 3 Gaussian quadrature points are used for each random variable, so m N = 3 2 simulations are required, where m is the number of Gaussian quadrature points and N is the number of random variables. Similarly, 9 simulations are required for NIPC2.

From the results of the calculations above, we conclude that the results of UT, SC and NIPC2 agree well. However the number of calculations required by UT increases fast with the number of random variables. The performances of SC and NIPC2 are better than UT.

Combined uncertainties

Table 3.11 presents the results of the combined uncertainties in the peak 10g average SAR value induced by x 1 and x 2 . The individual uncertainties are calculated by UT, SC and NIPC2, respectively. Since the two random variables are assumed to be independent, the correlation coefficient ρ 12 = 0, and the combined uncertainties are close to the values in Table 3.10, i.e., close to the results of two-dimensional UQ methods. Thus, when there are multiple random variables, combined uncertainty is an alternative and interesting approach to obtain uncertainty in the output. However, when the correlations between the variables are unknown, it may be possible to assume the maximum correlation, i.e., ρ 12 = 1, thereby leading to a relatively higher value for the uncertainties. For example, in this case, the increase is about 0.5% and 1% at 890 MHz and 1750 MHz, respectively. Table 3.12 shows the results of another example when there are three uncertain input parameters x 1 , x 2 and x 3 , where x 3 = 4.8 ± 10% is the relative permittivity of the Liquid Crystal Display (LCD) and also assumed to be a Gaussian distribution. The three random variables are assumed to be independent, i.e., the correlation coefficients are Uncertainty quantification using non-intrusive methods ρ 12 = 0, ρ 13 = 0 and ρ 23 = 0. The maximum values of uncertainties are given when the correlation coefficients are equal to 1.

Sensitivity analysis in SAR calculation

Since the non-intrusive PC using projection method (NIPC2) is efficient for UQ in SAR calculation, it is straightforward to calculate the sensitivity indices following the UQ procedure.

Conclusion

Since it is very hard to apply MCM for UQ and SA in SAR calculation, in this thesis, alternative computational efficient approaches are proposed for UQ and SA. Each method has its own limitations and drawbacks as mentioned above. For example, UT is not very efficient when there are multiple uncertain inputs in the calculation. The performance of the regression method based on the ordinary least-square minimization (NIPC1) is not good. The full truncation schemes in the regression approach have the drawback which is that they require a dramatically increasing number of model evaluations when the dimensionality (i.e., the number of input parameters in the model) increases. To circumvent this problem, sparse representations of the truncation have been studied in order to reduce the computational cost such as LARS that may provide a better solution [5]. NIPC2 and SC have better performances than the other approaches in our calculations. However, they face the curse of dimensionality problem when the number of uncertain inputs increases. Combining uncertainties which follows the clue of MoM is proposed to solve this problem. In the case studied herein, the random variables may be considered to be independent since the contribution for maximum correlation applying gPC expansion in 1D-FDTD simulation is then given. The results of the method are compared with MCM and NIPC.

Stochastic system formulation

Formulating stochastic systems aims at giving an existing deterministic model for an electromagnetic system in order to study the impact of uncertain input parameters to the system [START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF]. The deterministic governing equation can be formulated as follows.

Let Ω ⊂ R d , d = {1, 2, . . . } be a fixed physical domain, x = {x 1 , . . . , x d } be the coordinates. For the electromagnetic problem, the governing equation can be written as [START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF]:

(x, t)u(x, t) = R(x, t), (x, t) ∈ Ω × T (4.1)
where t ∈ [0, T ] is the time interval of interest, u(x, t) is the solution, and R(x, t) is the source term. generally represents a differential operator which involves differentiations in space/time. The solution can be u(x, t) = E(x, t) or u(x, t) = H(x, t) where E(x, t) and H(x, t) represent electric field and magnetic field, respectively. Setting boundary and initial conditions, the governing equation can be solved by conventional numerical methods, e.g., the FDTD method.

Generally, the computational methods for electromagnetic simulations are based on the fact that the input parameters such as permittivity ε, conductivity σ, etc., are assumed to be known precisely. When the input parameters are given, we can solve the governing equation (4.1) to obtain the solutions such as the electric field and magnetic fields.

However, in practice, some severe stochastic variations may exist in the input parameters, and the effect of these variations must be considered in the electromagnetic simulation. To quantify the effect of variations in input parameters is a meaningful but challenging task.

Let ξ = {ξ 1 , . . . , ξ N } ∈ R N represent independent random variables with the PDF W i (ξ), for i = 1, . . . , N . Then the joint PDF W (ξ) is given by:

W (ξ) = N i=1 W i (ξ i ), ∀ξ ∈ Γ (4.2)
where Γ represents the support of the random variables.

The uncertain input parameters must be correctly characterized in terms of random variables with appropriate distributions. It depends on the availability of detailed exper-imental data regarding these parameters. However, for certain parameters, experiments provide limited information about the variations in input parameters. When there is not sufficient information on the input parameters, the most straightforward way is to model the uncertain parameter as a normally or uniformly distributed random variable over the given range [START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF]. Following that, the governing equation (4.1) can be reformulated in a stochastic manner in the case of uncertain input parameters ξ in the electromagnetic system as follows:

(x, t, ξ)u(x, t, ξ) = R(x, t, ξ), (x, t, ξ) ∈ Ω × T × Γ (4.3)
where u(x, t, ξ) = H(x, t, ξ) or u(x, t, ξ) = E(x, t, ξ) are the solutions. H(x, t, ξ) is magnetic field and E(x, t, ξ) is electric field.

The Finite Difference Time Domain method

Considering that Maxwell's equations are given in differential or integral form, the CEM methods used to solve EMC problems can be divided into two categories: differential equation methods and integral equation methods. The differential equation methods include FDTD, TLM, FEM and Finite-Difference Frequency-Domain (FDFD) [START_REF] Joannopoulos | Photonic Crystals: Molding the Flow of Light[END_REF][START_REF] Christ | Three-dimensional finite-difference method for the analysis of microwave-device embedding[END_REF][START_REF] Albani | A numerical method based on the discretization of Maxwell equations in integral form[END_REF].

The integral equation methods include MoM and Fast Multiple Method (FMM) [START_REF] Walton | The Method of Moments in Electromagnetics[END_REF].

However, many CEM methods allow a solution to be calculated in either the frequency or the time domain. In this situation the CEM methods can be split into two groups: frequency domain and time domain methods. The first ones used in EMC include FDFD, FEM, MoM and FMM. The second ones include FDTD and TLM. The FDTD method is considered in this section.

In [7], Yee proposed a numerical method for solving electromagnetic problems in the time domain. This method, known as the FDTD method, is based on Maxwell's equations for electric and magnetic fields. The motive was to choose a geometry for spatially sampling the electric and magnetic field vector components which robustly represent both differential and integral forms of Maxwell's equations. The equations of 1D-FDTD are shown below.

The propagation of electromagnetic waves is determined by Maxwell's two curl equations. In their most general form these equations read as [21]:

▽ ×E = - ∂ B ∂ t (4.4) ▽ ×H = ∂ D ∂ t + J (4.5)
where E is the electric field, D the electric displacement field, H the magnetic field, B the magnetic flux density and J the electric current density. Considering a simple case where the electromagnetic fields are propagating through a linear, isotropic, homogeneous dielectric medium with a permittivity, permeability and conductivity denoted by ε, µ and σ, respectively. In such a medium the following relations hold:

D =
∂ H x ∂ t = 1 µ ∂ E y ∂ z - ∂ E z ∂ y (4.11) ∂ H y ∂ t = 1 µ ∂ E z ∂ x - ∂ E x ∂ z (4.12) ∂ H z ∂ t = 1 µ ∂ E x ∂ y - ∂ E y ∂ x (4.13)
and

∂ E x ∂ t = 1 ε ∂ H z ∂ y - ∂ H y ∂ z -σE x (4.14) ∂ E y ∂ t = 1 ε ∂ H x ∂ z - ∂ H z ∂ x -σE y (4.15) ∂ E z ∂ t = 1 ε ∂ H y ∂ x - ∂ H x ∂ y -σE z (4.16)
For a particular case, one considers an electromagnetic wave propagating in the negative x-direction, with its electric field component oriented in the z-direction, and its magnetic field component oriented in the y-direction. Then Maxwell's equations reduce to

∂ H y ∂ t = 1 µ ∂ E z ∂ x (4.17) ∂ E z ∂ t = 1 ε ∂ H y ∂ x -σE z (4.18)
The above two equations are hyperbolic partial differential equations. In [START_REF] Lax | Systems of conservation laws[END_REF], the authors proposed the use of a central derivative approximation to solve such equations. For a function f (x) the central derivative approximation takes the form

∂ f (x) ∂ x ≈ f (x + s/2) -f (x -s/2) s (4.19)
where s represents a small constant. This approximation is a finite difference approximation. Applying the finite difference approximation to equations (4.17 The interval along the x-axis can be split into points that are separated by ∆x, where the jth point is given by x j = j∆x. Similarly, time can be split into time steps ∆t, where the nth time step is denoted as t n = n∆t [21]. This forms a mesh of points in space and time. The fields E z (x j , t n ) are denoted as E n z ( j). The update equations are obtained 4.3 The Finite Difference Time Domain method 85 by applying this discretisation to the above equations, and they are

H n+1/2 y ( j + 1/2) = H n-1/2 y ( j + 1/2) + ∆t µ E n z ( j + 1) -E n z ( j -1) ∆x (4.26) E n+1 z ( j) = 1 -∆tσ/2ε 1 + ∆tσ/2ε E n z ( j) + ∆t ε(1 + ∆tσ/2ε) H n+1/2 y ( j + 1/2) -H n+1/2 y ( j -1/2) ∆x (4.27)
The update equations of FDTD are

H n+1/2 y ( j + 1/2) = H n-1/2 y ( j + 1/2) + γ E n z ( j + 1) -E n z ( j -1) (4.28) E n+1 z ( j) = αE n z ( j) + β H n+1/2 y ( j + 1/2) -H n+1/2 y ( j -1/2) (4.29)
where the parameters γ, α and β are given by:

γ = ∆t µ∆x (4.30) α = 1 -∆tσ/2ε
1 + ∆tσ/2ε (4.31)

β = ∆t ∆xε(1 + ∆tσ/2ε) (4.32)
The update equations (4.28) and (4.29) represent the FDTD algorithm for the propagation of electromagnetic waves along the x-axis. To use these equations the electric and magnetic field values need to be given initial values at time t = 0. Initially the field values are set at zero to represent free space.

To initialise an electromagnetic wave in the one-dimensional mesh a certain amount of electric field is added at certain points in the mesh over a number of time steps. The amount of field added at each time step, and the points at which the field is added depend on the input excitation that is being modelled. The input excitation for position j and time step n can be represented as E n s ( j). Inserting this input excitation into the FDTD update equation (4.29) yields

E n+1 z ( j) = αE n z ( j) + β(H n+1/2 y ( j + 1/2) -H n+1/2 y ( j -1/2)) + E n s ( j) (4.33)
The 1D-FDTD is relatively easy to implement into a computer code, and following the idea of 1D-FDTD, the equations of 3D-FDTD can be obtained. 3D-FDTD is summarized in Appendix A.

Generalized polynomial chaos in finite difference time domain simulation

The solutions of equation (4.3) E(x, t, ξ) and H(x, t, ξ) are regarded as second-order random processes, and they can be expanded by gPC as [21]:

H(x, t, ξ) = P i=0 h i (x, t)ψ i (ξ) (4.34)
σ 2 H = 〈H(x, t, ξ), H(x, t, ξ)〉 -〈H(x, t, ξ), 1〉 2 = P i=0 h 2 i 〈ψ 2 i 〉 -h 2 0 = P i=1 h 2 i 〈ψ 2 i 〉 (4.44) σ 2 E = 〈E(x, t, ξ), E(x, t, ξ)〉 -〈E(x, t, ξ), 1〉 2 = P i=0 e 2 i 〈ψ 2 i 〉 -e 2 0 = P i=1 e 2 i 〈ψ 2 i 〉 (4.45)
In 1D-FDTD update equations, the solutions H y and E z can be expanded onto polynomials ψ, and they are rewritten as

H n y ( j, ξ) = P i=0 h n i ( j)ψ i (ξ) (4.46) E n z ( j, ξ) = P i=0 e n i ( j)ψ i (ξ) (4.47)
where h n i and e n i are polynomial coefficients. To obtain the mean and uncertainty of the output fields, the field coefficients h n i and e n i must be calculated. The steps to calculate them are presented in the next section.

Considering that there are uncertain parameters in EMC simulation as mentioned above, the 1D-FDTD update equations are formed as follows [21]:

H n+1/2 y ( j + 1/2, ξ) = H n-1/2 y ( j + 1/2, ξ) + γ(E n z ( j + 1, ξ) -E n z ( j -1, ξ)) (4.48) E n+1 z ( j, ξ) = αE n z ( j, ξ) + β(H n+1/2 y ( j + 1/2, ξ) -H n+1/2 y ( j -1/2, ξ)) (4.49)
Substituting equations (4.46) and (4.47) into (4.48) and (4.49) yields the following equations [21]:

P i=0 h n+(1/2) i ( j + 1/2)ψ i (ξ) = P i=0 [h n-(1/2) i ( j + 1/2)ψ i (ξ) + γ(e n i ( j + 1) -e n i ( j -1))ψ i (ξ)] (4.50) P i=0 e n+1 i ( j)ψ i (ξ) = P i=0 [αe n i ( j)ψ i (ξ) + β(h n+1/2 i ( j + 1/2) -h n+1/2 i ( j -1/2))ψ i (ξ)] (4.51)
The electric field source E s can be put into equation (4.51), and the expansion of E s is:

E n s = P i=0 e n i,s ( j)ψ i (ξ) (4
( j + 1/2) = h n-1/2 k ( j + 1/2) + 1 〈ψ 2 k 〉 P i=0 [(e n i ( j + 1) -e n i ( j -1))〈γψ i ψ k 〉] (4.54) e n+1 k ( j) = e n+1 k,s ( j) + 1 〈ψ 2 k 〉 P i=0 [e n i ( j)〈αψ i ψ k 〉 + (h n+1/2 i ( j + 1/2) -h n+1/2 i ( j -1/2))〈βψ i ψ k 〉] (4.55)
If the material properties ε, µ and σ all depend on the uncertain parameter ξ, then γ, α and β are given by:

γ = ∆t µ(x, ξ)∆x (4.56) α = 1 -∆tσ(x, ξ)/2ε(x, ξ) 1 + ∆tσ(x, ξ)/2ε(x, ξ) (4.57) β = ∆t ∆xε(x, ξ)(1 + ∆tσ(x, ξ)/2ε(x, ξ)) (4.58)
The inner products 〈αψ i ψ k 〉, 〈βψ i ψ k 〉 and 〈γψ i ψ k 〉 must be calculated first. That can be done using numerical integration. Then the update equations are in a similar manner to the leapfrog scheme used in conventional FDTD [21].

The boundary condition is required to complete the scheme. In [21], an example of the scattering of electromagnetic fields in free space is discussed, the first-order Mur absorbing condition being chosen as a relevant boundary condition. The update Intrusive polynomial chaos method for uncertainty quantification equation for the electric field at the lower boundary reduces to

e n+1 k (0) = e n k (1) + ζ(1)(e n+1 k (1) -e n k (0)) (4.59)
The update equation at the upper boundary j = M is

e n+1 k (M ) = e n k (M -1) + ζ(M -1)(e n+1 k (M -1) -e n k (M )) (4.60)
where k = 0, . . . , M and ζ is defined as

ζ( j, ξ) = β( j, ξ)γ( j, ξ) -1 β( j, ξ)γ( j, ξ) + 1 (4.61)
When the coefficients are calculated by the above equations, the mean and variance of the magnetic field H n y ( j, ξ) are obtained as mentioned above.

Uncertainty quantification applied to a one dimensional example

One example of intrusive gPC in 1D-FDTD simulation with one uncertain input is given here. Fig. 4.1 shows the geometry of this example. In this 1D-FDTD example, a Gaussian pulse is introduced into the one-dimensional problem space. It contains a dielectric slab, which has an uncertain relative permittivity.

The relative permittivity of the slab follows a Normal distribution with mean ε r = 4 and standard deviation σ ε r = 0.4. The width of the slab b = 0.1m, and distance a is 0.3m. The output electric field E z is recorded at the centre of the problem space, which is E z (x = 0.5, t). The Gaussian pulse is propagated from the point x = 0.7m in this 1D problem space. The form of the input excitation is observation point is

E s (x = 0.7, t) = E 0 ex p - (t 0 -t) 2
E z (x = 0.5, t) = E 0 ex p - (x 0 /c 0 + t 0 -t) 2 2k 2 + R 12 ex p - (x 0 /c 0 + 2a/c 0 + t 0 -t) 2 2k 2 + ∞ n=1 T 12 R 2n-1 21 T 21 ex p - (x 0 + 2a)/c 0 + 2mb/c + t 0 -t 2 2k 2 (4.63)
where T 12 and R 12 represent the transmission and reflection coefficients for a wave passing from free space to the dielectric slab. T 21 and R 21 represent the transmission and reflection coefficients for a wave passing from the dielectric slab to free space. The speed of the electromagnetic wave in free space and the dielectric are represented by c 0 and c, respectively. The last term contains an infinite sum, which is approximated by taking the sum to n = 30.

The analytic time response and the time response obtained by FDTD simulation at the specified observation point x = 0.5m are presented by Fig. 4.2. The relative permittivity of the dielectric slab is given the mean value 4. The absorbing boundary condition used in FDTD is the first-order Mur absorbing boundary condition [21]. The mean and uncertainty of the output of this FDTD simulation are calculated by different UQ methods in the following section. The relative permittivity is related to ξ which follows a Normal distribution with mean equal to zero and standard deviation equal to The material parameter α, used in the FDTD Polynomial Chaos update equations (4.57) is

α =          1-∆tσ/2ε r ε 0 1+∆tσ/2ε r ε 0 0.1 ≤ x ≤ 0.2 1-∆tσ/2ε 0 1+∆tσ/2ε 0 other wise (4.65)
where ε 0 is the permittivity of free space, and in practice the conductivity σ is given a very low non-negative value. The material parameter β, used in the update equations (4.58) is

β =          ∆t ∆xε r ε 0 (1+∆tσ/2ε r ε 0 ) 0.1 ≤ x ≤ 0.2 ∆t ∆xε 0 (1+∆tσ/2ε 0 )
other wise (4.66)

Uncertainty quantification applied to a one dimensional example 93

When the chaotic expansion is terminated at P = 3, the Hermite polynomials required are

ψ 0 = 1 ψ 1 = ξ ψ 2 = ξ 2 -1 ψ 3 = ξ 3 -3ξ (4.67)
If a point in the space is inside the slab (0.1 ≤ x ≤ 0.2), the inner products 〈αψ i ψ k 〉 and 〈βψ i ψ k 〉 are calculated by integrals presented in Appendix (A.3). Numerical integrations are used to calculate the above integrals, and the method is given in [21].

If a point in the space is outside the slab (0.1 > x and 0.2 < x), the inner products 〈αψ i ψ k 〉 and 〈βψ i ψ k 〉 are calculated as

〈αψ i ψ k 〉 = α〈ψ i ψ k 〉 (4.68) 〈βψ i ψ k 〉 = β〈ψ i ψ k 〉 (4.69)
The material parameter γ, used in the FDTD Polynomial Chaos update equations (4.56) has no ξ dependence, and its inner products 〈γψ i ψ k 〉 reduce to

〈γψ i ψ k 〉 = γ〈ψ i ψ k 〉 (4.70)
The inner products calculated by these integrals are substituted into the update equations (4.54) and (4.55), which are used to calculate the coefficients of polynomial chaos. The infinite limits of the improper integrals are transformed to finite limits.

The mean and uncertainty of the electric field in time-domain are obtained as mentioned above, however, one is often more interested in the statistics of the signal in the frequency-domain. For this purpose, the electric field is transformed to the normalised electric field by taking the absolute value of this field in frequency-domain relative to the absolute value of the incident pulse in frequency-domain. It is so called normalised electric field as follows

E nor m ( f ) = {E z (t)} {E s (t)} = E z ( f ) E s ( f ) (4.71)
where {} denotes the Fourier transform, E z (t) and E z ( f ) are the electric field in the time-domain and frequency-domain, respectively, E s (t) and E s ( f ) are the incident pulse in time-domain and the one which is transformed to the frequency domain, respectively, and | | denotes the absolute value. Then, the polynomial coefficients e l nor m ( f ) of the normalised electric field, for l = 0, ..., P, are calculated as follows:

e l nor m ( f ) = P i=0 {e i (t)}ψ i (ξ) P i=0 {e i,s (t)}ψ i (ξ) , ψ l (ξ) (4.72)
where {e i (t)} and {e i,s (t)} are the Fourier transforms of the ith coefficient of the uncertain time-domain signal E z (t) and the source E s (t), respectively. Numerical integrations are required to obtain these polynomial coefficients, which means extra computational time.

Results and discussions

In the one-dimensional example above, the output of interest is the normalised electric field. The mean normalised electric field and uncertainty in the normalised electric field are calculated by MCM, NIPC and intrusive gPC, respectively. The results are presented in Fig. 4.3 and Fig. 4.4. The number of simulations for MCM to obtain the mean and uncertainty is 1000. For both NIPC and intrusive gPC expansion, the order of polynomial chaos is 3. For NIPC, the projection method is used. The numbers of simulations required for NIPC and intrusive gPC expansion are 3 and 1, respectively. Fig. 4.3 shows that the mean normalised electric field obtained by the three UQ methods agree well up to about 1000MHz. Comparing with the mean normalised electric field obtained by intrusive gPC expansion, the result predicted by NIPC is closer to the result obtained by MCM. However, at some particular frequencies, for example, around 2400MHz or 2700MHz, the mean values predicted by NIPC or intrusive gPC expansion become more different from the result obtained by MCM. To find out the reason, the relationship between the output (normalised electric field) and the input (relative permittivity) at different frequencies is observed, and the results are presented in Fig. 4.5. Since the uncertainty in the relative permittivity is ±10% which is not small, the relative permittivity changes from 3.6 to 4.4. At 120MHz and 500MHz, the relationship between the normalised electric field and the relative permittivity is linear. At 2400MHz, 2760MHz and 3000MHz, the relationship is non-linear. The result given by Fig. 4.5 indicates that if the relationship between the output and input is linear or in a relative linear fashion, the performance of PC expansion may have a better performance. In [21], the author indicates that small perturbations in the input parameters, in some of the CEM examples, caused frequency shifts in the frequency response of the output of interest. The shifts in the frequency can cause large changes in the resonant output of interest at certain frequencies. This means that the output of interest can have a highly non-linear dependence on the input parameters at particular frequencies.

The uncertainties quantified by MCM and NIPC are similar at the frequencies below 500MHz. At high frequencies, the uncertainty quantified by NIPC become more different from the result of MCM, which may be due to the non-linear relationship between the output and the input as observed in Fig. 4.5. In fact, when we give the uncertain input a smaller variance, e.g., ±3%, the uncertainties calculated by NIPC and MCM agree well up to about 2000MHz, which further explains that the relationship between the uncertain input and output influences the performance of PC expansion. Comparing with NIPC, the uncertainty of the output predicted by intrusive gPC expansion is more different from the results obtained by MCM.

The different ways to implement gPC expansion in NIPC and intrusive gPC may be the reason to explain the results presented in Fig. 4.3 and Fig. 4.4. For intrusive gPC expansion, at every update step which is similar to FDTD update step, the inner products are required, and the numerical integrations are performed. The error due to numerical integration may influence the result in every step. For NIPC, since it is a post-process method, only in the last step which forms the mean value and uncertainty of the output, the Gaussian quadrature is used. The error due to quadrature rule only exists in this last step. The results of every simulation required by NIPC is not influenced by it.

The advantage of intrusive gPC expansion is that it requires only one simulation to obtain the mean and uncertainty of the output, which is computationally efficient. For some complicated calculation such as SAR calculation, it is an important feature because one SAR calculation may take several hours. If the moderate estimates of the mean and uncertainty in EMC calculation are required, intrusive gPC expansion is available. Good absorbing boundary conditions may improve the performance of intrusive gPC expansion [21].

Conclusion

There are three errors related to gPC expansion (intrusive and non-intrusive gPC expansion): the first error due to the finite truncation of infinite expansion, the second error due to the quadrature rule, and the third error due to numerical integration. The first error exists in both intrusive gPC expansion and NIPC. Since the coefficients of PC expansion in NIPC are calculated by the projection method which uses the Gaussian quadrature, the error due to the quadrature rule exists in NIPC. The intrusive gPC expansion is related to the third error.

Intrusive gPC expansion requires fewer simulations than NIPC or MCM to obtain the mean and uncertainty of the output, which is an important advantage when the calculation takes much time. However, implementing gPC expansion in FDTD simulation is difficult, especially, for 3D-FDTD simulations with multiple uncertain input parameters in the calculation. When the governing equations of the problem take highly complicated forms, an explicit derivation of the new code may become impossible [START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF]. The single simulation run of intrusive gPC in FDTD can take long time for 3D-FDTD simulations. From several calculations' results, the author of [21] gives the conclusion that intrusive gPC expansion predicts uncertainties with a similar level of accuracy as the MoM, but it is computationally more expensive when the numerical integrations are implemented and cannot be used in certain situations. Until now, to our present knowledge, no idea is proposed for this method to provide a solution about how to quantify the uncertainty in the output of a simulation when there are multiple correlated uncertain input parameters.

Since SAR calculation is a complicated problem which needs to consider modifying the existing code and calculating correlations between input parameters, using the intrusive gPC expansion is a very difficult task. In this thesis, we do not provide results of quantifying the uncertainty in SAR calculation using this intrusive method.

Chapter 5

Conclusions and Future Work

Conclusions

Since it is very difficult to implement Monte Carlo Method (MCM) in Specific Absorption Rate (SAR) calculation due to the fact that MCM is computationally expensive, investigating computational efficient Uncertainty Quantification (UQ) methods is required. The main contribution of this thesis is investigating some computationally efficient non-intrusive UQ methods in SAR calculation. The main aim of the work is using acceptable time or number of simulations to evaluate the uncertainty in SAR calculation, which is similar as the result obtained by MCM.

The non-intrusive UQ methods applied are Unscented Transformation (UT), Stochastic Collocation (SC) and Non-Intrusive Polynomial Chaos (NIPC), and their implementations are as straightforward as MCM, requiring only repetitive executing of a deterministic code at a set of predefined points in the stochastic space. These methods are satisfactory under specific conditions. UT gives good results for uncertainty analysis, yet the number of simulations it requires increases fast with the dimension of random vectors. It is not appropriate to implement UT when there are multiple uncertain input parameters in the calculation. When the collection points of SC and the polynomial coefficients of NIPC are calculated by quadrature rule, there are errors existing in the results of the two methods due to this quadrature rule. Using the regression method to evaluate the polynomial coefficients in NIPC, implementing oversampling is a good choice to improve the performance of the method, however, the efficiency of the method is decreased. The error due to the finite truncation of infinite expansion exists in SC and NIPC.

To overcome the curse of dimensionality in these UQ methods, combining uncertainties is proposed. It is an alternative approach to evaluate the overall uncertainty in the case of multiple uncertain input parameters in the simulation. The drawback of the method is that the variation in the input parameters should be small enough, otherwise, this method cannot be accurate.

The algorithm of applying gPC expansion in FDTD is presented in chapter 4, and it was investigated by other researchers and presented in their publications. This method has demonstrated excellent computation efficiency compared to the MCM. It can provide the statistical information of the outputs by a single simulation run, however, NIPC method and other non-intrusive methods require running the simulation several times to obtain the uncertainty and mean of the output.

One disadvantage of intrusive gPC expansion is that it requires a reformulation of the deterministic numerical solver for handling inner product integrals, which is difficult to implement. For example, the governing equations of SAR calculation take highly complicated forms and an explicit derivation of the new code for SAR calculation become difficult. In addition to the difficulty in the derivation of the new code, using numerical integrations to estimate the inner products makes this method suffering from large memory requirements, and increases its computational time. The single simulation run required by intrusive gPC method can take long time for a 3D-FDTD simulation.

Future Work

In fact, the more time-consuming the calculation is, the more obvious the advantage of intrusive gPC expansion becomes because comparing with other UQ methods, it requires a lesser number of calculations to obtain the uncertainty in the output. Since SAR calculation may take much time, future investigations targeted at applying intrusive gPC expansion in SAR calculation should be carried out following the work of this thesis.

For this challenging task, the first step could be applying the intrusive gPC expansion to SAR calculation with one uncertain input in the calculation, which requires rewriting the FDTD code for SAR calculation. In the second step, the situation of multiple uncertain inputs in SAR calculation should be considered, these uncertain inputs being correlated or uncorrelated. The problem of correlated uncertain inputs in intrusive gPC expansion has not been investigated.

The work in this thesis concentrated on quantifying the parametric uncertainties in the output of SAR simulations, but the model form uncertainty is not considered.

For three-dimensional FDTD (3D-FDTD), the Maxwell's equation in rectangular coordinates are as follows [START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF]:

µ ∂ H x ∂ t + σ m H x = ∂ E y ∂ z - ∂ E z ∂ y (A.1) µ ∂ H y ∂ t + σ m H y = ∂ E z ∂ x - ∂ E x ∂ z (A.2) µ ∂ H z ∂ t + σ m H z = ∂ E x ∂ y - ∂ E y ∂ x (A.3) µ ∂ E x ∂ t + σ e E x = - ∂ H y ∂ z + ∂ H z ∂ y (A.4) µ ∂ E y ∂ t + σ e E y = - ∂ H z ∂ x + ∂ H x ∂ z (A.5) µ ∂ E z ∂ t + σ e E z = - ∂ H x ∂ y + ∂ H y ∂ x (A.6)
where µ is permeability, ε is permittivity and σ e and σ m are the electric and the magnetic conductivities, respectively. x, y and z denote the coordinates.

Intrusive gPC expansion in 3D-FDTD Let (i, j, k) = (i∆x, j∆ y, k∆z). The update equations of FDTD become: The inner products in these equations must be evaluated firstly, and it is a pre-process.

H
Then the updating equations can be solved in a way similar to the traditional FDTD method.

3. 7

 7 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, NIPC1 without oversampling and NIPC1 with two times oversampling and three times oversampling, respectively, at 900MHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, NIPC1 without oversampling and NIPC1 with two times oversampling and three times oversampling, respectively, at 1800MHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 (a) CAD-based model of the commercially available mobile phone, (b) the mobile phone model placed against the right cheek position of the SAM phantom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 SAR distributions of the numerical phone model given in Fig. 3.9, and the line as presented is the one considered with Fig. 3.12. (a) SAR distribution at 890MHz, (b) SAR distribution at 1750MHz. . . . . . . . . 3.11 Measurement results of the SAR distributions of the mobile phone Neo FreeRunner. (a) SAR distribution at 890MHz, (b) SAR distribution at 1750MHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvii 3.12 SAR values in the middle line of SAR distributions given in Fig. 3.10. (a) SAR values in the middle line of the SAR distribution of 890 MHz, (b) SAR values in the middle line of the SAR distribution of 1750 MHz. . . . 3.13 Sensitivity analysis by PC based Sobol' indices in SAR calculation, and the order of PC r=2, (a) Sensitivity indices at 890MHz, (b) Sensitivity indices at 1750MHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 A one dimensional problem space containing a dielectric slab. A Gaussian pulse is excited at ,the resulting fields are observed at ×. . . . . . . . . 4.2 The analytic time response and the time response obtained by FDTD simulation of the one dimension example (The relative permittivity of the dielectric slab is given the mean value 4). . . . . . . . . . . . . . . . . . 4.3 Mean normalised electric field calculated using the three different UQ methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Uncertainty in the normalised electric field calculated using the three different UQ methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Normalised electric field with different relative permittivity at different frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33 Fig. 2 . 2

 3322 Fig. 2.1 Dual-band mobile phone model.

Fig. 2 . 3 Fig. 2 . 4

 2324 Fig. 2.3 In the case of only random variable x 1 in the calculation: the mean and standard deviation of the output calculated by MCM and second-order UT, respectively, for 900MHz

2 Fig. 2 . 5 2 Fig. 2 . 6

 225226 Fig. 2.5 In the case of two random variables : the mean and standard deviation of the output calculated by MCM and second-order UT, respectively, for 900MHz

Fig. 2 . 7

 27 Fig. 2.7 In the case of one random variable, the mean and standard deviation of the output calculated by MCM, second-order UT and fourth-order UT, respectively, for 900MHz

  It has been introduced by Choi et al. (2004b) [63], Berveiller (2005) [64] based on early results in ocean engineering by Tatang (1995) [65], Tatang et al. (1997)

  .35) By integrating equation (3.29), the total variance can thus be decomposed into partial variances:

  (a) The mean of output, and the random variable is x 1 .(b) The standard variance of output, and the random variable is x 1 .

Fig. 3 .

 3 Fig. 3.1 With one random variable in the calculation: the mean and standard deviation of the output calculated by MCM, second-order UT, SC and NIPC2 respectively, at 900MHz

Fig. 3 . 2

 32 Fig. 3.2 With one random variable in the calculation: the mean and standard deviation calculated the output by MCM, second-order UT, SC and NIPC2 respectively, at 1800MHz

Fig. 3 . 4

 34 Fig. 3.4 Three times oversampling is implemented in NIPC1: the mean and standard deviation of the output at 900MHz are calculated by MCM and NIPC1, respectively, and NIPC1 uses three different sets of sampling points represented as NIPC1 Case1, NIPC1 Case2, and NIPC1 Case3, respectively.

  (a) The mean of output, and the random variables are x 1 and x 2 . (b) The standard deviation of output, and the random variables are x 1 and x 2 .

Fig. 3 . 5

 35 Fig. 3.5 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, second-order UT, SC and NIPC2 respectively, at 900MHz.

  (a) The mean of output, and the random variables are x 1 and x 2 . (b) The standard deviation of output, and the random variables are x 1 and x 2 .

Fig. 3 . 6

 36 Fig. 3.6 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, second-order UT, SC and NIPC2 respectively, at 1800MHz.

  (a) The mean of output, and the random variables are x 1 and x 2 . (b) The standard deviation of output, and the random variables are x 1 and x 2 .

Fig. 3 . 7

 37 Fig. 3.7 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, NIPC1 without oversampling and NIPC1 with two times oversampling and three times oversampling, respectively, at 900MHz.

Fig. 3 . 8

 38 Fig. 3.8 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, NIPC1 without oversampling and NIPC1 with two times oversampling and three times oversampling, respectively, at 1800MHz.

  3 provides the dielectric properties of the tissue equivalent liquid for each frequency band. The SAR distributions at 890 MHz and 1750 MHz of the numerical model are shown in Fig. 3.10 (a) and (b), respectively. The SAR distributions obtained by experiments or measurements undertaken by the laboratory leading in [10] are given in Fig. 3.11.

Fig. 3 .

 3 Fig. 3.9 (a) CAD-based model of the commercially available mobile phone, (b) the mobile phone model placed against the right cheek position of the SAM phantom.

Fig. 3 .

 3 Fig. 3.10 SAR distributions of the numerical phone model given in Fig. 3.9, and the line as presented is the one considered with Fig. 3.12. (a) SAR distribution at 890MHz, (b) SAR distribution at 1750MHz.

Fig. 3 .

 3 Fig. 3.11 Measurement results of the SAR distributions of the mobile phone Neo FreeRunner. (a) SAR distribution at 890MHz, (b) SAR distribution at 1750MHz.

Fig. 3 .

 3 Fig. 3.12 SAR values in the middle line of SAR distributions given in Fig. 3.10. (a) SAR values in the middle line of the SAR distribution of 890 MHz, (b) SAR values in the middle line of the SAR distribution of 1750 MHz.

Fig. 3 .

 3 13 shows the sensitivity indices of Sobol decomposition, and Sobol' indices are calculated by PC expansion with two uncertain inputs. The two random variables are x 1 and x 2 as mentioned above. The overall indices are equal to 100%.The order of PC is 2. SU 1 and SU 2 are the sensitivity indices of x 1 and x 2 , respectively, and SU 12 is the sensitivity index of the combination of x 1 and x 2 . The sensitivity indices show that the relative permittivity of the antenna support x 1 has a dominant impact on the calculation results. When the number of uncertain input parameters increases, it becomes more difficult to calculate Sobol' indices by PC expansion.

Fig. 3 .

 3 Fig. 3.13 Sensitivity analysis by PC based Sobol' indices in SAR calculation, and the order of PC r=2, (a) Sensitivity indices at 890MHz, (b) Sensitivity indices at 1750MHz.

  )

Fig. 4 .

 4 Fig. 4.1 A one dimensional problem space containing a dielectric slab. A Gaussian pulse is excited at ,the resulting fields are observed at ×.

Fig. 4 . 2

 42 Fig. 4.2 The analytic time response and the time response obtained by FDTD simulation of the one dimension example (The relative permittivity of the dielectric slab is given the mean value 4).

Fig. 4 . 3

 43 Fig. 4.3 Mean normalised electric field calculated using the three different UQ methods.

Fig. 4 . 4

 44 Fig. 4.4 Uncertainty in the normalised electric field calculated using the three different UQ methods.

Fig. 4 . 5

 45 Fig. 4.5 Normalised electric field with different relative permittivity at different frequencies.

A. 3 〈αψ 0 ψ 2 〉1〈βψ 0 ψ 2 〉〈βψ 0 ψ 3 〉 = 〈βψ 3 ψ 0 2 - 3 〈βψ 2 ψ 3 〉 = 〈βψ 3 ψ 2

 322302332 The inner products of intrusive polynomial chaos expansion〈αψ 0 ψ 1 〉 = 〈αψ 1 ψ 0 〉 = 〈αψ 2 ψ 0 〉 -∆tσ/2(4 + 0.4ξ)ε 0 1 + ∆tσ/2(4 + 0.4ξ)ε 0 (ξ 3 -3ξ) ex p(-4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) 〈βψ 0 ψ 1 〉 = 〈βψ 1 ψ 0 〉 4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) = 〈βψ 2 ψ 0 〉 4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) (ξ 2 -1) ex p(-4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) 4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) ξ(ξ 2 -1) ex p(-4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) ξ(ξ 3 -3ξ) ex p(-4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) (ξThe inner products of intrusive polynomial chaos expansion115 4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) (ξ 2 -1)(ξ 3 -3ξ)ex p(-4 + 0.4ξ)ε 0 )(1 + ∆tσ/2(4 + 0.4ξ)ε 0 ) (ξ 3 -3ξ) 2 ex p(-

2 Uncertainty quantification using the Unscented Transformation (UT

  Aims and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Outline of thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 The basic idea of UT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 The formulation of two-dimensional UT . . . . . . . . . . . . . . . . 2.3.3 Higher-order UT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Example 1: using second-order UT for UQ with one random variable in the calculation . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Example 2: using second-order UT for UQ with two random variables in the calculation . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Example 3: using fourth-order UT for UQ with one random variable in the calculation . . . . . . . . . . . . . . . . . . . . . . . . . .

	xxii	Table of contents
	2.3.1

Introduction 1.1 Uncertainty in SAR calculation . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Specific absorption rate (SAR) . . . . . . . . . . . . . . . . . . . . . 1.1.2 Finite Integration Technique (FIT) and SAR calculation in CST . 3 1.1.3 Uncertainty in modelling . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Probability approach for uncertain quantification . . . . . . . . . . . . . . 8 1.2.1 Random variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.2 Random vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Overview of the methods related to uncertainty quantification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Method of moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.2 Quadrature method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4 ) 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Random sampling and stratified sampling . . . . . . . . . . . . . . 20 2.2.2 Latin hypercube sampling (LHS) . . . . . . . . . . . . . . . . . . . . 2.3 Unscented transformation (UT) . . . . . . . . . . . . . . . . . . . . . . . . .

3 Uncertainty quantification using non-intrusive methods 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Stochastic collocation (SC) method . . . . . . . . . . . . . . . . . . . . . . . 3.3 Polynomial chaos (PC) expansion . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Homogeneous chaos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 The generalized polynomial chaos (gPC) method . . . . . . . . . . 3.3.3 Using non-intrusive polynomial chaos for UQ . . . . . . . . . . . . 3.3.4 Sensitivity analysis (SA) using PC expansion . . . . . . . . . . . . . 3.4 Combined uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Uncertainty quantification using a simple mobile phone model . 3.5.2 Numerical model for SAR calculation . . . . . . . . . . . . . . . . . 3.5.3 Uncertainty in SAR calculation . . . . . . . . . . . . . . . . . . . . . 3.5.4 Combined uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.5 Sensitivity analysis in SAR calculation . . . . . . . . . . . . . . . . . 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Intrusive polynomial chaos method for uncertainty quantification 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Stochastic system formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 The Finite Difference Time Domain method . . . . . . . . . . . . . . . . . . 4.4 Generalized polynomial chaos in finite difference time domain simulation 4.5 Uncertainty quantification applied to a one dimensional example . . . . 4.6 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table of contents xxiii 5 Conclusions and Future Work

 of 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Finite difference time domain in three dimensions . . . . . . . . . . . . . . 103 A.2 Generalized polynomial chaos in three-dimension finite-difference timedomain simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 A.3 The inner products of intrusive polynomial chaos expansion . . . . . . . . 111 Synoptic view of modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The framework of uncertainty analysis . . . . . . . . . . . . . . . . . . . . . 2.1 Dual-band mobile phone model. . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The return losses of the dual-band mobile phone as considered. . . . . . With one random variable in the calculation: the mean and standard deviation calculated the output by MCM, second-order UT, SC and NIPC2 respectively, at 1800MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Appendix A Intrusive gPC expansion in 3D-FDTD A.1 Appendix B Glossary List of figures 3.2 3.3 Non oversampling is implemented in NIPC1: the mean and standard List of figures deviation of the output at 900MHz are calculated by MCM and NIPC1, respectively, and NIPC1 uses three different sets of sampling points represented as NIPC1 Case1, NIPC1 Case2, and NIPC1 Case3, respectively. 3.4 Three times oversampling is implemented in NIPC1: the mean and standard deviation of the output at 900MHz are calculated by MCM 1.1 xxvi and NIPC1, respectively, and NIPC1 uses three different sets of sampling
	points represented as NIPC1 Case1, NIPC1 Case2, and NIPC1 Case3,
	References

2.3 In the case of only random variable x 1 in the calculation: the mean and standard deviation of the output calculated by MCM and second-order UT, respectively, for 900MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 In the case of only random variable x 1 in the calculation: the mean and standard deviation of the output calculated by MCM and second-order UT, respectively, for 1800MHz . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 In the case of two random variables : the mean and standard deviation of the output calculated by MCM and second-order UT, respectively, for 900MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 In the case of two random variables : the mean and standard deviation of the output calculated by MCM method and second-order UT, respectively, for 1800MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 In the case of one random variable, the mean and standard deviation of the output calculated by MCM, second-order UT and fourth-order UT, respectively, for 900MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 In the case of one random variable, the mean and standard deviation of the output calculated by MCM, second-order UT and fourth-order UT, respectively, for 1800MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 With one random variable in the calculation: the mean and standard deviation of the output calculated by MCM, second-order UT, SC and NIPC2 respectively, at 900MHz . . . . . . . . . . . . . . . . . . . . . . . . . . respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, second-order UT, SC and NIPC2 respectively, at 900MHz. . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 With two random variables in the calculation: the mean and standard deviation of the output calculated by MCM, second-order UT, SC and NIPC2 respectively, at 1800MHz. . . . . . . . . . . . . . . . . . . . . . . . .

Table 2 .

 2 1 The results for the case of x 2 in the calculation.

	Random	UQ method	Mean and uncertainty Mean and uncertainty
	variable		at 900 MHz	at 1800 MHz
	x 2	MCM	0.895±1.06%	1.718±1.35%
		Second-order UT	0.895±1.06%	1.721±1.34%

Table 3 .

 3 1 Correspondence of the orthogonal polynomials and random variables for different Askey-chaos.

		Random variables x Orthogonal polynomials ψ i	Support
	Continuous	Normal	Hermite	(-∞, +∞)
		Gamma	Laguerre	[0,∞)
		Beta	Jacobi	[a,b]
		Uniform	Legendre	[a,b]
	Discrete	Poisson	Charlier	{0,1,2,. . . }
		Binomial	Krawtchouk	{0,1,. . . , m}
		Negative Binomial	Meixner	{0,1,2,. . . }
		Hypergeometric	Hahn	{0,1,. . . , m}

  and SU 1,2,...,k are the sensitivity indices corresponding to a set of random variables {x 1 , x 2 , ...,x k } T ∈ R k . σ 2

y

is the variance of the output given in equation (3.23), ψ s and s are presented as in equation

(3.19)

, however, at least one s i in the set

Table 3 .

 3 2 Dielectric properties of different components of the phone

	Component	Relative permittivity	Conductivity
	Antenna	1.00	PEC
	Antenna support	2.33	0.01
	Battery	1.00	PEC
	Battery connectors	1.00	PEC
	Casing	3.00	0.01
	LCD	4.80	0.01
	PCB	1.00	PEC
	Receiver	1.00	PEC
	Speaker	1.00	PEC
	Speaker connectors	1.00	PEC
	Vibrator	1.00	PEC
	PEC: perfect electrical LCD: liquid crystal PCB: printed circuit
	conductor	display	board
	Table 3.3 Dielectric properties of tissue equivalent liquid
	Frequency	890MHz 1750MHz
	Relative permittivity	41.50	40.00
	Conductivity [S/m]	0.97	1.40
	Mass density [kg/m 3 ]	1000	1000

Table 3 .

 3 4 SAR calculation at 890 MHz using different normalizations.

		x 1	SAR10g SAR1g
			[W/kg] [W/kg]
	SAR normalized with	1.93 1.0542 1.4545
	respect to same accepted 2.33 1.0500 1.4467
	power (250mW): SAR acc	2.73 1.0191 1.4034
	SAR normalized with	1.93 0.8538 1.1777
	respect to same stimulated 2.33 0.7148 0.9848
	power (250mW): SAR st i	2.73 0.5683 0.7826
	Table 3.5 SAR calculation at 1750 MHz using different normalizations.
		x 1	SAR10g SAR1g
			[W/kg] [W/kg]
	SAR normalized with	1.93 0.5787 0.9866
	respect to same accepted 2.33 0.6336 1.0800
	power (125mW): SAR acc	2.73 0.6512 1.1070
	SAR normalized with	1.93 0.3968 0.6765
	respect to same stimulated 2.33 0.5906 1.0067
	power (125mW): SAR st i	2.73 0.6258 1.0638
	Table 3.6 Uncertainty of SAR st i at 890MHz
	u SAR st i	Uncertainty of Uncertainty of
	SAR10g[W/kg] SAR1g[W/kg]
	Quantified by SC	16.34%	16.41%
	Combined uncertainty: u c	15.98%	15.48%
	Table 3.7 Uncertainty of SAR st i at 1750MHz
	u SAR st i	Uncertainty of Uncertainty of
	SAR10g[W/kg] SAR1g[W/kg]
	Quantified by SC	18.45%	18.38%
	Combined uncertainty: u c	12.02%	12.00%

Table 3 .

 3 8 UQ in SAR calculation by one dimensional UQ methods,

	Random	Method Number	SAR10g	SAR10g
	Variable		of Sim-	[W/kg]	[W/kg]
			ulations	890(MHz)	1750(MHz)
		UT	3	1.046±1.13% 0.605±4.58%
	x 1 =	SC	3	1.044±1.60% 0.627±5.02%
	2.33±10% NIPC2	3	1.041±1.77% 0.625±4.82%
		UT	3	1.049±0.25% 0.633±1.11%
	x 2 =	SC	3	1.050±0.34% 0.633±1.47%
	3.00±10% NIPC2	3	1.049±0.35% 0.633±1.61%

Table 3 .

 3 9 UQ in SAR calculation by NIPC1.

	Random Number	SAR10g	SAR10g
	Variable of Sim-	[W/kg]	[W/kg]
		ulations	890(MHz)	1750(MHz)
		6	1.047±3.12% 0.640±13.05%
		9	1.041±2.53% 0.639±10.16%
	x 1	12	1.043±2.25% 0.636±8.75%
		15	1.044±2.23% 0.628±5.27%
		18	1.043±2.22% 0.628±6.30%

Table 3 .

 3 10 UQ in SAR calculation by two dimensional UQ methods.

	Random Method Number	SAR10g	SAR10g
	Variables		of Sim-	[W/kg]	[W/kg]
			ulations	890(MHz)	1750(MHz)
		UT	36	1.046±1.09% 0.615 ±4.65%
	x 1 , x 2	SC	9	1.044±1.58% 0.627±4.93%
		NIPC2	9	1.040±1.79% 0.624±5.54%

Table 3 .

 3 11 Combined uncertainty: u c

	Correlation Method Number	SAR10g	SAR10g
	Coefficient		of Sim-	[W/kg]	[W/kg]
			ulations 890(MHz) 1750(MHz)
		UT	6	1.22%	4.62%
	ρ 12 = 0	SC	6	1.63%	5.27%
		NIPC2	6	1.80%	5.02%
		UT	6	1.38%	5.61%
	ρ 12 = 1	SC	6	1.94%	6.47%
		NIPC2	6	2.11%	6.41%

Table 3 .

 3 12 Combined uncertainty: u c

	Correlation Method Number	SAR10g	SAR10g
	Coefficients:		of Sim-	[W/kg]	[W/kg]
			ulations 890(MHz) 1750(MHz)
	ρ 12 = ρ 13 =	UT	9	1.16%	4.57%
	ρ 23 = 0	SC	9	1.63%	5.19%
		NIPC2	9	1.88%	5.03%
	ρ 12 = ρ 13 =	UT	9	1.21%	5.63%
	ρ 23 = 1	SC	9	1.97%	6.55%
		NIPC2	9	2.20%	6.49%

2 Generalized polynomial chaos in three-dimension finite-difference time-domain simulation

  Generalized polynomial chaos in three-dimension finite-difference time-domain simulation 107The method mentioned in chapter 4 can be extended to the 3D-FDTD simulation. The solutions E and H of the 3D-FDTD can be expanded by polynomials, and they are given by[START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF]:

	n+ 1 2 x A.1 Finite difference time domain in three dimensions i, j + 1 2 , k + 1 2 = A hxh H n x i, j + 1 2 , k + 1 2 + A hx e y E n y i, j + 1 2 , k + 1 -E n y i, j + + A hx ez E n z i, j + 1, k + 1 2 -E n y i, j, k + 1 2 where A hxh = H n+ 1 2 z i + 1 2 , j + 1 2 , k = A hzh H n z i + 1 2 , j + 1 2 , k + A hze y E n y i + 1, j + 1 2 , k -E n y i, j + 1 2 + A hze x E n x i + 1 2 , j + 1, k -E n x i + 1 2 , j, k , k 1 2 , k where A hzh = 2µ -σ m ∆t 2µ + σ m ∆t Intrusive gPC expansion in 3D-FDTD (A.7) 105 (A.9) E n+1 y i, j + 1 2 , k = A e y e E n y i, j + 1 2 , k + A e yhx H n+ 1 2 x i + 1 2 , j, k + 1 2 -H n+ 1 2 x i + 1 2 , j, k -1 2 + A e yhz H n+ 1 2 z i + 1 2 , j + 1 2 , k -H n+ 1 2 z i -1 2 , j + 1 2 , k (A.11) where 2ε -σ e ∆t P Intrusive gPC expansion in 3D-FDTD e n+1 x m (i + 1 2 , j, k) = P i=0     e n x i (i + 1 2 , j, k)〈A e x e ψ i (ξ)ψ m (ξ)〉 +(h n+ zi (i + 1 2 , j + 1 2 , k) -h n+ 1 2 zi (i + 1 2 , j -1 2 , k))〈A e xhz ψ i (ξ)ψ m (ξ)〉 +(h n+ y i (i + 1 2 , j, k + 1 2 ) -h n+ 1 2 y i (i + 1 2 , j, k -1 2 ))〈A e xh y ψ i (ξ)ψ m (ξ)〉     (A.27) A.2 E n x (i, j, k, ξ) = e n (A.13) x i (i, j, k)ψ i (ξ) A e y e = i=0 2ε + σ e ∆t P 2µ -σ m ∆t 2µ + σ m ∆t A hx e y = A hze y = -2∆t (2µ + σ m ∆t)∆x A e yhx = 2∆t E n y (i, j, k, ξ) = i=0 e n (A.14) y i (i, j, k)ψ i (ξ) (2ε + σ e ∆t)∆z 2∆t (2µ + σ m ∆t)∆z 2∆t A hze x = 2∆t (2µ + σ m ∆t)∆ y P 2∆t A e yhz = -(2ε + σ e ∆t)∆x E n z (i, j, k, ξ) = i=0 e n zi (i, j, k)ψ i (ξ) (A.15) e n+1 y m (i, j + 1 , k) 2 A hx ez = -(2µ + σ m ∆t)∆ y H n+ 1 2 y i + 1 2 , j, k + 1 2 = A h yh H n y i + 1 2 , j, k + 1 2 + A h y e x E n x i + 1 2 , j, k + 1 -E n x i + 1 2 E n+1 x i + 1 2 , j, k = A e x e E n x i + 1 2 , j, k + A e xhz H n+ 1 2 z i + 1 2 , j + 1 2 , k -H n+ 1 2 z i + 1 2 , j -1 , k E n+1 z i, j, k + 1 2 = A eze E n z i, j, k + 1 2 H n x (i, j, k, ξ) = P i=0 h n x i (i, j, k)ψ i (ξ) (A.16) = P i=0     e n y i (i, j + 1  2 , k)〈A e y e ψ i (ξ)ψ m (ξ)〉 +(h n+ x i (i + 1 2 , j, k + 1 2 ) -h n+ 1 2 x i (i + 1 2 , j, k -1 2 ))〈A e yhx ψ i (ξ)ψ m (ξ)〉 +(h n+ zi (i + 1 2 , j + 1 2 , k) -h n+ 1 2 zi (i -1 2 , j + 1    (A.28) 2 , k))〈A e yhz ψ i (ξ)ψ m (ξ)〉 P 2 + A e xh y H n+ 1 2 y i + 1 2 , j, k + 1 2 -H n+ 1 2 y i + 1 2 , j, k -1 (A.10) + A ezh y H n+ 1 2 y i + 1 2 , j, k + 1 2 -H n+ 1 2 y i -1 2 , j, k + 1 H n y (i, j, k, ξ) = i=0 h n (A.17) y i (i, j, k)ψ i (ξ) 2 2 P , j, k + A h y ez E n z i + 1, j, k + 1 2 -E n z i + 1, j, k + 1 2 (A.8) where A h yh = 2µ -σ m ∆t 2µ + σ m ∆t A h y e x = -2∆t (2µ + σ m ∆t)∆z A h y ez = 2∆t (2µ + σ m ∆t)∆x where + A ezhx H n+ 1 2 x i, j + 1 2 , k + 1 2 -H n+ 1 2 x i, j -1 2 , k + H n z (i, j, k, ξ) = h n zi (i, j, k)ψ i (ξ) (A.18) 1 i=0 (A.12) 2 2ε -σ e ∆t A e x e = 2ε + σ e ∆t A e xhz = where Substituting equations (A.13)-(A.18) back into equations (A.7)-(A.12) to obtain : 2∆t (2ε + σ e ∆t)∆ y A e xh y = -2∆t (2ε + σ e ∆t)∆z A eze = 2ε -σ e ∆t 2ε + σ e ∆t A ezh y = 2∆t (2ε + σ m ∆t)∆x A ezhx = -2∆t (2ε + σ e ∆t)∆ y P i=0 h n+ 1 2 x i (i, j + 1 2 , k + 1 )ψ i (ξ) e n+1 zm (i, j, k + 1 ) 2 2 = P i=0        A hxh h n x i (i, j + 1 2 , k + 1 2 )ψ i (ξ) +A hx e y (e n y i (i, j + 1 2 , k + 1) -e n y i (i, j + 1 2 , k))ψ i (ξ)        = P i=0     e n zi (i, j, k + 1  2 )〈A eze ψ i (ξ)ψ m (ξ)〉 +(h n+ y i (i + 1 2 , j, k + 1 2 ) -h n+ 1 2 y i (i -1 2 , j, k + 1 2 ))〈A ezh y ψ i (ξ)ψ m (ξ) +(h n+ x i (i, j + 1 2 , k + 1 2 ) -h n+ 1 2 x i (i, j -1 2 , k + 1 2 )A ezhx ψ i (ξ)ψ m (ξ)〉    (A.29)

A.

+A hx ez (e n zi (i,

j + 1, k + 1 2 ) -e n y i (i, j, k + 1 2 ))ψ i (ξ)

(a) The mean of output calculated by MCM and NIPC1 which is without oversampling. (b) The standard deviation of output calculated by MCM and NIPC1 which is without oversampling.

(a) The mean of output calculated by MCM and NIPC1 which is with three times oversampling. (b) The standard deviation of output calculated by MCM and which is with three times oversampling.

(a) The mean of output, and the random variables are x 1 and x 2 . (b) The standard deviation of output, and the random variables are x 1 and x 2 .

Acknowledgements

Chapter 4

Intrusive polynomial chaos method for uncertainty quantification

Introduction

Before NIPC be proposed and studied, gPC expansion is used intrusively for UQ at first in the field of CEM. Since the gPC expansion is applied into a simulation code (such as FDTD), the methods are so-called intrusive gPC.

In [START_REF] Chauvière | Computational modelling of uncertainty in time-domain electromagnetics[END_REF] work involving the implementation of the gPC expansion into a higherorder discontinuous Galerkin solution of Maxwell's equations has been proposed. The gPC method can accurately quantify the uncertainty in the output, while being more computationally efficient than the MCM. This work, however, only estimated the output uncertainty due to a single uncertain input parameter. Then, in [21], the authors proposed a new algorithm to apply the gPC expansion. They applied it in FDTD simulation in order to quantify the uncertainty in CEM simulation output. The update equations of the new algorithm follow a form similar to the one of conventional FDTD. The results of the output uncertainty which is induced by three independent uncertain input parameters evaluated by the new method for a simple CEM simulation example are then given. Otherwise, Shen [START_REF] Shen | Efficient numerical techniques for uncertainty analysis in electromagnetic simulations[END_REF] applied the new algorithm for various types of CEM simulation examples such as backscattered field and radar-cross-section in scattering applications to analyse the impact of uncertainties on simulation results, and the method was proved efficient.

In the present chapter, the new algorithm of applying the gPC expansion in FDTD simulation is introduced in detail. Before introducing the method, basics such as stochastic system formulation and FDTD method are summarised. An example of

101

The model form uncertainty is the uncertainty in the output of a simulation due to the uncertainty in how well the model used represents reality, which is the most important type of uncertainty. Much work should be carried out to determine the model form uncertainty in SAR simulations.

Appendix A

Intrusive gPC expansion in 3D-FDTD A.1 Finite difference time domain in three dimensions

Intrusive gPC expansion in 3D-FDTD

A.2 Generalized polynomial chaos in three-dimension finite-difference time-domain simulation 109

Applying the inner product to both sides of equation (A. 19)-(A.24) with polynomial ψ m and 0 ≤ m ≤ P yields

Intrusive gPC expansion in 3D-FDTD

2 ) 2π dξ (A.34)

2 ) 2π dξ (A.37)