Higher Order Modelling for Computational Electromagnetics - - PowerPoint PPT Presentation

Higher Order Modelling for Computational Electromagnetics Computational Electromagnetics

Roberto D. Graglia Dipartimento di Elettronica Politecnico di Torino – Italy

e-mail: roberto.graglia@polito.it

1

WHERE IS TORINO (TURIN)?

• In the heart of Europe

e ea

• u ope
• The taste of industry
• The industry of taste
• The shape of technology
• Resources for development

TORINO

Turin, the “Capital of the Alps,” is well known as the home of the Shroud of Turin, the headquarters

• f

automobile manufacturers Fiat Lancia and Alfa Romeo manufacturers Fiat, Lancia and Alfa Romeo (Torino is “the Automobile Capital

• f

Italy”), and as host of the 2006 Winter

• Olympics. It is one of Italy's main industrial

centers, being part of the famous “industrial triangle,” along with Milan and Genoa. Torino was founded by the Romans in the first century BC, it was the capital of the Duchy Torino was founded by the Romans in the first century BC, it was the capital of the Duchy

• f Savoy from 1563, then of the Kingdom of Sardinia ruled by the Royal House of Savoy

and finally the first capital of unified Italy (1861). Th it tl h t f It l ' b t i iti ll d i h th The city currently hosts some of Italy's best universities, colleges, academies, such as the Polytechnic University of Turin. Prestigious and important museums, such as the Egyptian Museum (the 2nd in the world after the Cairo one) and the Mole Antonelliana are also found in the city. Turin's several monuments and sights make it one of the world's top y g p 250 tourist destinations, and the tenth most visited city in Italy in 2008.

Francesco Vercelli (22/10/1883 - 24/11/1952) Geofisico piemontese, laureato a Torino nel 1908 in fisica e nel 1909 in matematica, si occupò principalmente del mare e delle altre acque Amedeo Avogadro (09/08/1776 09/07/1856) Scienziato piemontese che studiò le proprietà dei gas con Amedeo Avogadro (09/08/1776 - 09/07/1856) Scienziato piemontese che studiò le proprietà dei gas con l'enunciazione del principio che porta il suo nome Gustavo Colonnetti (08/11/1886 - 20/03/1968) Il fondatore dell’Istituto di Metrologia del CNR di Torino Bernardino Drovetti (1776 - 1852) Archeologo e diplomatico piemontese "padre" del Museo Egizio di Torino Cesare Lombroso (06/11/1835 - 19/10/1909) Psichiatra, antropologo e criminologo del XIX secolo Giuseppe Luigi Lagrange (25/01/1736 - 10/04/1813) Matematico torinese famoso in tutta l’Europa di fine Giuseppe Luigi Lagrange (25/01/1736 - 10/04/1813) Matematico torinese famoso in tutta l Europa di fine Settecento Giovanni Antonio Amedeo Plana (06/11/1781 - 20/01/1864) Astronomo, fu il vero fondatore dell'Osservatorio Astronomico di Torino. La sua fama è specialmente legata agli scritti su la "Teoria della Luna" Galileo Ferraris (30/10/1847 - 07/02/1897) Elettrotecnico, studiò i motori a correnti alternate e fondò a Torino la prima scuola italiana per ingegneri elettrotecnici Giovanni Schiaparelli (14/03/1835 - 04/07/1910) Uno dei più prolifici astronomi italiani del XIX secolo, famoso p ( ) p p , per le sue osservazioni della superficie di Marte Primo Levi (31/07/1919 - 11/04/1987) Scrittore ebreo piemontese celebre per aver narrato l’orrore dei campi di concentramento, chimico di professione, la sua opera è un interessante esempio di connubio fra letteratura e i scienza Rita Levi Montalcini (22/04/1909) Scienziata neurobiologa piemontese vincitrice nel 1986 del premio Nobel per la medicina

The Politecnico has 29.300 students studying on 96 courses (22 Bachelor's degree courses; 31 Master of Science courses; 23 Doctorates and 20 specialization courses). 18 of them are ; p ) held in English. In the academic year 2010/2011 the Politecnico had around 4.800 students in the first year; in 2009 over 4.000 students graduated with a Master of Science or a Bachelor's Degree. Each year, between lectures, laboratories and practical exercises there are 170 000 hours of teaching There is a staff of over 900 lecturers and researchers and are 170.000 hours of teaching. There is a staff of over 900 lecturers and researchers, and around 875 administration staff. There are 5 Schools, 1 Graduate School, 18 Departments and 7 Interdepartmental Centers. The income in the 2010 forecast balance is 380 million Euros (in 1990 the figure was 52 million). The Ministero dell'Istruzione Università e Ricerca or M.I.U.R. (Ministry for Universities Education and Research) contributes around 129 million Euros.

Higher Order Modelling for Computational Electromagnetics Computational Electromagnetics

Roberto D. Graglia Dipartimento di Elettronica Politecnico di Torino – Italy

e-mail: roberto.graglia@polito.it

6

SUMMARRY Boundary value problems & numerical approaches 8 - 13 20 minutes Differential and integral operators 14 - 24 Conclusion of the first part 25 - 27 Higher order modelling 28 - 83 35 minutes Results for high-order vector-bases 84 - 91 Results for high-order singular vector-bases 92 - 104

7

C id d ti & bl Considered equations & problems

• We deal with Electromagnetic problems modelled by

We deal with Electromagnetic problems modelled by Maxwell’s equations.

• The techniques described have applications to many

practical problems, e.g.: EMC & EMP; Shielding radiation from printed circuits; Microwave hazards; Electromagnetic radiation from and penetration into vehicles, aircrafts, radiation from and penetration into vehicles, aircrafts, ships, missiles; Antennas near ground; Design

• f

frequency selective surfaces for reflector antennas and radomes; Radar scattering; Propagation in

• ptical

radomes; Radar scattering; Propagation in

• ptical

components & systems; etc.

8

9

Numerical approach to boundary value problems Numerical approach to boundary value problems B d l bl Boundary value problem Integral formulation

Mi ed or h brid form lation

Differential formulation

Mixed or hybrid formulation

Numerical methods to represent operator equations in terms of matrix equations

10

Numerical methods (for linear problems) Numerical methods (for linear problems)

• in 3 slides -

For a linear, inhomogeneous equation: Introduce a set of testing functions {wm}, and a suitable inner product

11

Then obtain a matrix equation for a square, non-singular matrix with The solution may be exact The solution

• r approximate depending on the choice of {fn} and {wm}

if f G l ki ’ th d may be exact

12

if wn=fn → Galerkin’s method

How to choose {f } {w }? How to choose {fn}, {wm}?

• The set {fn} (linearly independent functions) has to

approximate f reasonably well;

• The set {fn} could be chosen in order to satisfy the

BCs (FEM applications) ; BCs (FEM applications) ;

• The set {wm} (linearly independent functions) has to

be chosen so that <w g> depends on relatively be chosen so that <wm, g> depends on relatively independent properties of g;

• The choice depends on the desired accuracy of the

solution;

• The ease of evaluation of matrix elements;
• The realization of a well-conditioned matrix [ℓm,n].

13

Differential & Integral operators Differential & Integral operators

(a) boundary conditions; (b) well posedness

Example: uniform transmission line [α=xb=π/2]

• differential operator

p The boundary conditions

• integral operator

MUST be enforced g p Boundary conditions are automatically SATISFIED

14

This reads as an integral equation of the 1st kind

Thi diff i l bl i ll d h l i f

(b) well posedness

This differential problem is well posed = the solution f is “weakly” affected by small variations of the source g

i t f DIFFERENTIAL bl t t 1 input of DIFFERENTIAL problem input perturbed by adding a gaussian-pulse 0.1 in height 0 25 0.3

• utput

perturbed output 0.8 0 15 0.2 0.25 0.4 0.6 g(x) 0.05 0.1 0.15 f(x) 0.2

• 0.05

15

0.2 0.4 x/ 0.2 0.4

• 0.1

x/

(b) well posedness

This integral problem is not well posed = the solution g is “badly” affected by small variations of the source f

1 input of INTEGRAL problem input perturbed by adding a gaussian-pulse 0.001 in height 7

• utput

perturbed output 0.8 1 5 6 0 4 0.6 f(x) 4 5 g(x) 0.2 0.4 2 3

16

0.2 0.4 x/ 0.2 0.4 1 x/

Why on Earth is this integral problem not well-posed? not well posed?

17

Because the kernel of this integral problem is too smooth → it is even continuous at x=x’ !!!!!!!

(well-posedness is guaranteed whenever the operator is not compact)

0.2 0.4 x

 )

• 0.2

l Kernel K(x, 0 6

• 0.4

Integral K(x, x'=/4) K(x, x'=3/8) K(x, x'=/2) 18

• 0.5

0.5 1

• 0.6

x/ envelope for the maxima of K(x,x')

In case of infinite domain (i e

• ∞ <x< ∞) the

In case of infinite domain (i.e., ∞ <x< ∞), the Green function relative to this differential problem with “natural” boundary conditions is f ti f ( ’) a function of (x-x’)

K=[cos(α|x-x’|)-sin(α|x-x’|)= [ ( | |) ( | |) =√2 Re{exp(-j π/4) exp(-j α|x-x’|)}

Convolution type operators permit one to reduce by appropriate use

• f space Fourier transforms the number of integrations required by

the used numerical technique → No time to discuss this here.

19

Integral equations for EM applications Integral equations for EM applications, and their kernels

∫ ζ g(r)- ∫K(r,r’) g(r’) dr’=f(r) f(r) is known and g(r) is unknown ( ) g( ) ζ=0 → eq. of the first kind (EFIE for PEC); ζ≠0 → eq of the second kind (HFIE) ζ≠0 → eq. of the second kind (HFIE).

The kernel must be singular to get well-posed first kind eqs.

Very Luckily, in EM, one usually has singular kernels, at least of logarithmic kind for 2D bl lik 1/ f 3D bl

20

problems or like 1/r for 3D problems

Representations of Fields in Terms of Potentials (Frequency-Domain)

• Solutions of Maxwell’s equations can be expressed in terms of potentials:

21

Potentials, Lorenz Gauges, and Continuity Equations

• In homogeneous, isotropic, unbounded regions
• Vector and scalar potentials are related by the Lorenz gauge:
• Charge densities are related to current densities by continuity equations:

22

Representations of Fields in Terms of Potentials (Frequency-Domain)

• Solutions of Maxwell’s equations can be expressed in terms of potentials:

23

Uniqueness of Solutions to Maxwell’s Equations

• Fields satisfying

in unbounded regions with specified on boundaries and radiation conditions at infinity are unique

• Fields satisfying Maxwell’s eqs. in bounded regions with lossy media

h i R ( / β) ≠ 0 i specified on boundaries and radiation conditions at infinity are unique.

• r having Re(α / β) ≠ 0 are unique.
• Uniqueness fails at certain discrete frequencies in lossless regions having

lossless surface impedance boundary conditions [Re(α / β)=0]; these frequencies

24

• ss ess su ace

peda ce bou da y co d t o s [ e(α / β) 0]; t ese eque c es correspond to the cavity resonant frequencies of the bounded region.

Generalization of the previous transparencies yields to the following

• To invert the differential operator is to solve the problem

→ to do this one has to know (find) the Green function

yields to the following

→ to do this one has to know (find) the Green function G of the given problem;

• If G is known, then the solution is simply obtained by

, p y y integration;

• Unfortunately, for all problems of practical interest G is

unknown → numerical approaches to differential/integral problems are required because of the ignorance of the Green function; ignorance of the Green function;

• The numerical problem could be well- or bad-posed; this

depends on the properties of the operator; p p p p ;

25

• In general, in integral formulations, use of the proper

( k ) G i id d b ti t l (unknown) G is avoided by resorting to general principles (e.g., equivalence principle), or to integral theorems (e.g., Green theorems); ( g , );

• Integral formulations guarantee the BCs and all the

physical properties of the solutions, as a matter of fact, in many cases, the proper integral equations are

• btained only at the moment of imposing the BCs (in

the finite region); the finite region);

• This is why IE are widely used in EM (since radiation

conditions / far-field conditions are automatically y satisfied);

• Because of the Green kernel, the matrix approximation
• f the problem is given by a dense matrix.

26

FEM MOM

Differential formulation Integral formulation

☺

Non linear problems can be dealt with Applications to linear problems are well known BCs have to be enforced BCs automatically satisfied

☺

Infinite domains: there is some problems Infinite domains: no problem

☺ ☺ Sparse matrices

Dense matrices (problems in dealing with very large problems)

☺

The integrals to be evaluated are The singularities of the kernel d i l i t ti it

☺

g simple render numerical integration quite difficult Complex geometry are easily dealt with Higher order models have been introduced more recently with introduced more recently

27

Higher order modelling Higher order modelling

Required: Required:

• To better represent the geometry of the problem.
• To reduce the L2 (least squares norm) of the

l i i f i solution error on regions of interest.

• It usually improves the convergence of the results and/or reduces the

size of the numerical matrices.

Higher order modelling

 Regarding the geometry of the problem; d t i t th t f th t need to approximate the curvature of the geometry; need to approximate small details in multiscale problems.

• use subsectional bases thereby subdividing the geometry into

y g g y simple (curved) elements of small size.

29

Use subsectional bases thereby subdividing the geometry into simple (curved) elements of small size into simple (curved) elements of small size.

30

Pablo Picasso: "Retrato de Ambroise Vollard” (1910).

Regarding the geometry of the problem

Contribution from specular points (in the optical limit).

Points of reflection on the body at which the angle of incidence is l t th l f fl ti l ti t th b ti i t equal to the angle of reflection relative to the observation point.

e g the RCS σ of a metallic sphere of radius a is σ=πa2

31

e.g., the RCS σ of a metallic sphere of radius a is σ=πa

In the high-frequency range, the PEC-sphere “models” shown here poorly model the commonly used “sphere benchmark” (see here poorly model the commonly used sphere-benchmark (see below) NASA Communications Satellite E h P j t (1960 1969) Echo Project (1960-1969)

32

Regarding the geometry of the problem

Travelling waves (contribution from points with high curvature). e.g., estimation of the near-nose-on cross-sections of long, thin bodies.

33

Regarding the geometry of the problem

Creeping waves

• f importance for the analysis in the shadow region

(the importance of the creeping waves depends on the dimension in wavelength of the body).

34

dimension in wavelength of the body).

Higher order modelling

 Regarding the e pansion/testing f nctions sed in the  Regarding the expansion/testing functions used in the numerical application. reduce the L2 (least squares norm) of the solution error on regions of interest; i f th lt improve convergence of the results; reduce the size of the numerical matrices.

• use vector functions of high polynomial order.

35

Results from a 1995 paper (16 years old): nose-on incidence on open circular wg. inner length=30 λ, inner diameter=5 λ, wall thickness 0.5 λ, profile length=66.5 λ Magnitude (LHS) and phase in radians (RHS) of the current functions for a truncated circular waveguide at incidence along the z-axis, obtained by solving a system with 816 unknowns.

36

z axis, obtained by solving a system with 816 unknowns.

How to g How to get cur t curved ed (i.e (i.e., distor ., distorted) elements ted) elements. How to g How to get cur t curved ed (i.e (i.e., distor ., distorted) elements ted) elements.

37

Source: Zinkiewicz’s book

How to g How to get cur t curved / ed / distor distorted elements ted elements.

• Any curved cell is obtained by mapping a “parent cell” into the
• bject domain.
• The parent cell for triangular patches is a rectilinear triangle.
• In object space, the 3 edges are the zero-coordinate lines ξ1, ξ2,

ξ ξ3=0.

• ξi is the area coordinate ξi = Ai /AT
• Dependency relations: ξ1 + ξ2 + ξ3 = 1

Dependency relations: ξ1 + ξ2 + ξ3 1

How to g How to get cur t curved / ed / distor distorted elements ed elements.

• The rectilinear triangle on the right-hand side is

r = ξ1 r1+ ξ2 r2+ ξ3 r3 ξ1 1 ξ2 2 ξ3 3

• Where ξ1, ξ2, ξ3 are 3 linear “shape functions,” each associated

to a different corner node of the cell.

Cur Curved elements; 2 ed elements; 2nd

nd or

• rder distor

der distortion ion

• Shape function: 2 ξ3 (ξ3 – ½) for node

Shape function: 2 ξ3 (ξ3 ½) for node

• Shape function: 4 ξ1 ξ3 for node
• The shape functions of the other nodes are obtained by cyclic

permutation of the subscripts permutation of the subscripts.

• Interpolatory polynomials easily provide the shape

functions.

Cur Curved elements; 2 ed elements; 2nd

nd or

• rder distor

der distortion ion

Shape functions: Mapping: Mapping:

Cur Curved elements; 2 ed elements; 2nd

nd or

• rder distor

der distortion ion

Shape functions:

Dependency relation, unitary basis vectors, and Jacobian and Jacobian

Mapping: Mapping:

Interpolatory polynomials for shape functions for shape functions. Use Lagrange interpolation polynomials written in g g p p y terms of interpolatory polynomials of Silvester.

Interpolatory polynomials of Silvester

• Use polynomials of degree i in ξ, where ξ is in the interval

[0 1] [0, 1].

• The parameter p indicates the number of uniform subintervals

p p into which the interval is divided.

• The polynomial is unity at ξ=i/p and has zeros at ξ=0 1/p 2/p
• The polynomial is unity at ξ=i/p and has zeros at ξ=0, 1/p, 2/p,

…, (i-1)/p

  



1 ) ( 1

1

k

i

          

 

1 1 ) ( ! ) , ( i p i k p i p R

i k i

  

Silvester polynomials: example for p=1

• There are 2 polynomials: one of 0th and one of 1st degree.

Pl d h l i l h i l [0 1]

• Please draw the polynomials on the interval [0, 1].

  ) , 1 ( ; 1 ) , 1 (

1

R R    1 ) , 1 ( ; 1 ) , 1 (

1

R R            



 

1 1 ) ( ! 1 ) , (

1

i p i k p i p R

i k i

     0 1 i

Silvester polynomials: example for p=1

• There are 2 polynomials: one of 0th and one of 1st degree.

Pl d h l i l h i l [0 1]

• Please draw the polynomials on the interval [0, 1].

     ) , 1 ( ; 1 ) , 1 (

1

R R

Silvester polynomials: example for p=2

• Three polynomials: one of 0th, one of 1st and one of 2nd degree.

Pl d h l i l h i l [0 1]

• Please draw these polynomials on the interval [0, 1].

    ) 1 2 ( 2 ) 2 ( ; 2 ) 2 ( ; 1 ) 2 ( R R R            1 2 1 ) , 2 ( ; 1 ) , 2 ( ; 1 ) , 2 (

2 1

R R R             



 

1 1 ) ( ! 1 ) , (

1

i p i k p i p R

i k i

 

Silvester polynomials: example for p=2

• Please draw these 3 polynomials on the interval [0, 1].

2 1 ) 1 2 ( 2 ) , 2 ( ; 1 2 ) , 2 ( ; 1 ) , 2 (

2 1

           R R R

Silvester polynomials: example for p=3

• 4 polynomials;
• ne of 0th; one of 1st; one of 2nd, one of 3rd degree.

Pl d th l i l th i t l [0 1]

• Please draw these polynomials on the interval [0, 1].

; 3 ) 3 ( ; 1 ) 3 ( R R           ) 2 3 )( 1 3 ( 3 ) 3 ( ; ) 1 3 ( 3 ) 3 ( ; 1 ) , 3 ( ; 1 ) , 3 (

3 2 1

R R R R               3 2 1 ) , 3 ( ; 2 1 ) , 3 (

3 2

R R            



 

1 1 ) ( ! 1 ) , (

1

i p i k p i p R

i k i

    1 i

Silvester polynomials: example for p=3

• Please draw these 4 polynomials on the interval [0, 1].

; 3 ) 3 ( ; 1 ) 3 (      R R 3 2 1 ) 2 3 )( 1 3 ( 3 ) , 3 ( ; 2 1 ) 1 3 ( 3 ) , 3 ( ; 1 ) , 3 ( ; 1 ) , 3 (

3 2 1

                R R R R 3 2 1 ) ( 2 1 ) (

3 2

    

Scalar Lagrangian interpolation on the g g p canonical elements: the segment.

• There are two ξ coordinates: ξ

ξ

• There are two ξ coordinates: ξ1, ξ2 .
• These are dependent coordinates (ξ1+ξ2 =1)

αij(ξ1, ξ2 )=Ri(p, ξ1)Rj(p, ξ2 ) with i+j=p is a p-th order Lagrangian polynomial interpolating points within a segment whose normalized coordinates (ξ1, ξ2 ) are (i/p, j/p).

51

Interpolation on a segment, example for p=1 p g , p p

• There are two interpolatory polynomials of 1st degree

There are two interpolatory polynomials of 1 degree.

• αij(ξ1, ξ2 )=Ri(p, ξ1)Rj(p, ξ2 ), (with i+j=1).
• Please draw the polynomials on the interval ξ1=[0 1]
• Please draw the polynomials on the interval ξ1 [0, 1].

) 1 ( ) 1 ( ) (         R R

2 2 1 1 2 1 01 1 2 1 1 2 1 10

) , 1 ( ) , 1 ( ) , ( ) , 1 ( ) , 1 ( ) , (                 R R R R

2 2 1 1 2 1 01

) , ( ) , ( ) , (     

52

Interpolation on a segment, example for p=1 p g , p p

) 1 ( ) 1 ( ) (         R R

2 2 1 1 2 1 01 1 2 1 1 2 1 10

) , 1 ( ) , 1 ( ) , ( ) , 1 ( ) , 1 ( ) , (                 R R R R

53

Interpolation on a segment, example for p=2 p g , p p

• Example for p=2: there are three interpolatrory

polynomials of 2nd degree polynomials of 2nd degree.

• αij(ξ1, ξ2 )=Ri(p, ξ1)Rj(p, ξ2 ), (with i+j=2).
• Please draw the polynomials on the interval ξ1=[0, 1].

Please draw the polynomials on the interval ξ1 [0, 1].

) 1 2 ( 2 ) , 2 ( ) , 2 ( ) , (

1 1 2 1 2 2 1 20

          R R 2 2 ) , 2 ( ) , 2 ( ) , ( 2 ) , 2 ( ) , 2 ( ) , (

2 1 2 1 1 1 2 1 11 2 1 2 2 1 20

              R R R R 2 ) 1 2 ( 2 ) , 2 ( ) , 2 ( ) , (

2 2 2 2 1 2 1 02

          R R

54

Interpolation on a segment for p=2 p g p

2 ) 1 2 ( 2 ) , 2 ( ) , 2 ( ) , (

1 1 2 1 2 2 1 20

          R R ) 1 2 ( 2 ) 2 ( ) 2 ( ) ( 2 2 ) , 2 ( ) , 2 ( ) , (

2 2 2 1 2 1 1 1 2 1 11

                   R R R R 2 ) , 2 ( ) , 2 ( ) , (

2 2 1 2 1 02

       R R

55

Scalar Lagrangian interpolation

• n a segment
• n a segment.
• For the p-th order Lagrangian polynomial

interpolation of a segment there are (p+1) interpolating polynomials of p-th degree.

56

Parameteriza meterization of tion of the points on a triang the points on a triangle le.

Recall that the coordinates of a rectilinear triangle may be parameterized as

,

3 3 2 2 1 1

r r r r      

p where ri is a vertex position vector for vertex i.

,

3 3 2 2 1 1

r r r r     

57

Parametriza metrization of ion of the points on a triang the points on a triangle le.

A Lagrange parametrization of order q for a curvilinear triangle can be expressed in terms of interpolating polynomials of the Silvester form as:

q k j i q R q R q R

q

     ) ( ) ( ) (    r r

where a triple indexing scheme is used to label the position

q k j i q R q R q R

k j k j i i ijk

    



), , ( ) , ( ) , (

3 2 , , 1

   r r

where a triple indexing scheme is used to label the position vector rijk interpolating the point with normalized coordinates (ξ1, ξ2, ξ3)=(i/q, j/q, k/q). (ξ1, ξ2, ξ3) ( q, j q, q)

Geometry description f Geometry description for a r a triang triangle le. Geometry description f Geometry description for a r a triang triangle le.

Normalized coordinates related by ξ1 + ξ2+ ξ3 =1 Ed d i d f h “i d d ” di ξ d ξ Edge vectors derived from the “independent” coordinates ξ1 and ξ2 ℓi=∂r/∂ ξi , i=1,2, from which the edge vectors that follow are found: from which the edge vectors that follow are found:

Edg Edge v vecto ectors s of

• f a triang

a triangle le. ℓ 1 = - ℓ 2, ℓ ℓ 1 ℓ 2 = ℓ 1, ℓ 3 = ℓ 2 - ℓ 1

Gr Gradient v adient vector ectors of s of a triang a triangle le. The gradient vectors are determined from the edge The gradient vectors are determined from the edge vectors as: “ξ = (n x ℓ ) /J “ξi = (n x ℓ i ) /J h ℓ 1 ℓ 2/J i th it t l t th where n = ℓ 1 x ℓ 2/J is the unit vector normal to the triangle while J= |ℓ 1 x ℓ 2| is the Jacobian.

Edg Edge and g and gradien adient v vecto ectors of s of a trian a triangle le.

The above definitions for edge, (height) and gradient vectors apply for triangles on curved surfaces with an extended interpretation Triangle tangent to a curvilinear triangle at a point. The curvilinear and (rectilinear) tangent triangles have the same element coordinates, Jacobian, edge vectors, and height vectors at the point of tangency.

Why to use cur Why to use curved ( ed (i.e i.e., distor , distorted) elements ted) elements? Why to use cur Why to use curved ( ed (i.e i.e., distor , distorted) elements ted) elements?

• the same model (that is, the same database) used by

structural /mechanical engineers can be used → this is a structural /mechanical engineers can be used → this is a big plus!!!

• better approximations of boundaries are obtained

with very small error in volume/surface/line models (no rescaling)

• usually a smaller number of unknowns is needed whenever

usually a smaller number of unknowns is needed whenever higher-order expansion functions are used → shorter computation time

63

Why to use cur Why to use curved ( ed (i.e i.e., distor , distorted) elements ted) elements? Why to use cur Why to use curved ( ed (i.e i.e., distor , distorted) elements ted) elements?

• mixed approaches (FEM + MoM) are facilitated because

the expansion functions of the two methods can be chosen the expansion functions of the two methods can be chosen within the “same” set.

• there is the possibility to construct singular functions able

to model singular current or field behaviors (though applications of this are rather new). But, is there any problem?

64

But, is there any problem? But, is there any problem?

• In FEM applications you now have to use quadrature to

evaluate the matrix entries, whereas for non-distorted cells there are often closed form results possible.

• For MoM applications people are already using quadrature
• For MoM applications people are already using quadrature,

but the treatment of the Green's function singularity might be different with curved cells.

• The main “cost” of curved cells is

(1) the additional data structure necessary in the model, (2) h ddi i l i i i f h d (2) the additional computation arising from the quadrature.

65

Vector functions

• We first consider the lowest order functions.
• Then we move on and consider higher-order vector

functions.

66

The family of generalized triangle basis functions (for current representation)

Source: Don Wilton’s notes

67

Div Div-Conf Conf functions on a triangular element functions on a triangular element Div Div Conf Conf functions on a triangular element functions on a triangular element

Divergence-conforming functions of the Nedelec type maintain only normal continuity across element maintain only normal continuity across element Boundaries (they do not prescribe tangential continuity).

68

Div Div-Conf Conf functions on a triangular element functions on a triangular element Div Div Conf Conf functions on a triangular element functions on a triangular element

They eliminates spurious solution of the EFIE operator while discarding the highest order degrees of freedom while discarding the highest order degrees of freedom associated with the nullspace of the divergence operator in the EFIE operator.

69

Div Divergence-conf ence-conforming bases on rming bases on triang triangles les.

1 



1 , 1

2 3 3 2 1

      J 

• Zeroth-order bases.
• Three vector basis

functions of first order

 ,

1

3 1 1 3 2

      J 

functions of first order.

• They have constant

normal and linear tangential (CN/LT)

 .

1

1 2 2 1 3

      J 

tangential (CN/LT) components at element edges.

J

70

Curl-conf l-conforming ba

• rming bases on tr

ses on triang iangles les.

: ˆ

 

    n

• Zeroth-order bases

Th t b i f ti

,

2 3 3 2 1

        

• Three vector basis functions
• f first order.
• They have constant

, ,

3 1 1 3 2 2 3 3 2 1

               

tangential and linear normal (CT/LN) components at element edges

.

1 2 2 1 3

        

(prove this by use of eq. (48) of 1997 paper).

71

Compl Completenes teness of

• f z

zeroth or th order triangula der triangular bases r bases (w (we consider only the cur e consider only the curl-conf conforming ones).

• rming ones).

(w (we consider only the cur e consider only the curl conf conforming ones).

• rming ones).

,

2 3 3 2 1

         , ,

3 1 1 3 2 2 3 3 2 1

             .

1 2 2 1 3

        

• The basis set is incomplete to first order since 6 degrees of

The basis set is incomplete to first order since 6 degrees of freedom are required to model linear variations in two independent vector components on a surface.

72

Compl Completenes teness of

• f z

zeroth or th order triangula der triangular bases r bases (w (we consider only the cur e consider only the curl-conf conforming ones).

• rming ones).

(w (we consider only the cur e consider only the curl conf conforming ones).

• rming ones).

,

2 3 3 2 1

         , ,

3 1 1 3 2 2 3 3 2 1

             .

1 2 2 1 3

        

• To render the bases first-order complete one must include the

p curl-free combinations:

,

1 1 

  , ,

2 2 1 1

          

73

.

1 2 2 1

      

Compl Completenes teness of

• f z

zeroth or th order triangula der triangular bases r bases (w (we consider only the cur e consider only the curl-conf conforming ones).

• rming ones).

(w (we consider only the cur e consider only the curl conf conforming ones).

• rming ones).

,

2 3 3 2 1

         , ,

3 1 1 3 2 2 3 3 2 1

             .

1 2 2 1 3

        

• Completeness to zeroth-order is proved by noticing that the

following linear combinations are able to represent two independent basis vectors on a 2D element (verify this by yourself and express “ξ3).

,

1 3 2

        

74

.

2 2 3

     

Compl Completenes teness of

• f z

zeroth or th order triangula der triangular bases r bases (w (we consider only the cur e consider only the curl-conf conforming ones).

• rming ones).

(w (we consider only the cur e consider only the curl conf conforming ones).

• rming ones).

,

2 3 3 2 1

         , ,

3 1 1 3 2 2 3 3 2 1

             .

1 2 2 1 3

        

The Nedelec conditions also require completeness of the The Nedelec conditions also require completeness of the curl to the same order as the bases. Completeness of the curl to zeroth-order follows from (verify this by use of (50, 53)

• f the 1997 paper)
• f the 1997 paper).

3 2 1 ˆ 2       n

75

3 , 2 , 1 ,      



n J

Higher or Higher order v der vector bases ector bases Higher or Higher order v der vector bases ector bases

76

Higher or Higher order bases der bases Higher or Higher order bases der bases

• One can construct higher order bases complete to order p by

f i h d f h d b i h l forming the product of zeroth-order bases with complete polynomial factors of order p.

• The set of polynomials factors used may take one of several

different forms chosen for convenience.

j i p s r :

• us

inhomogene p t s r : s homogeneou

s j r i t s r

       , , ,

3 2 1

     p k j i : al hierarchic p t s r :

• ry

interpolat        , ) , ) ) )

3 2 1 3 2 1

      , , ( H (p, R (p, R (p, R

ijk t s r

77

I t l t t f ti Interpolatory vector functions

• For interpolatory vector functions the key idea is to use shifted

p y y polynomials of Silvester to move interpolation points away from two of the edges-those along with the tangential (for curl- conf.) and the normal (for div.-conf.) components of the 0th- ) ( ) p

• rder basis factor vanish.

Interpolation nodes for curl- or divergence-conforming bases on triangular

• elements. Only nodes in basis subset Ω1

i j k or Λ1 i j k for p = 3 are shown.

78

Interpolation nodes for curl- or divergence-conforming bases on quadrilateral elements. Only nodes in basis subset Ω3

ik; jℓ or Λ3 ik; jℓ for

p = 2 are shown.

79

p

How to move the interpolation points i id h ll inside the cell

• We use the shifted polynomial of Silvester

i l i l f d ( ) i ξ h ξ i i h

) 1 , ( ) , ( ˆ

1

p p R p R

i i

 



 

• It is a polynomial of degree (i-1) in ξ, where ξ is in the

interval [0, 1] again subdivided in to p uniform subintervals.

• The polynomial is unity at ξ=i/p and has zeros at ξ=1/p,

2/p, …, (i-1)/p.

• There is no zero at ξ=0

There is no zero at ξ 0.

• For i¥1 one has

) ( ˆ ) (    p R p p R  ) , ( ) , (   p R i p R

i i



How to move the interpolation points How to move the interpolation points inside the cell

• Shifted polynomial of Silvester

1 2 ) ( 1

1

i k

i

    





 1 1 1 1 2 ) ( )! 1 ( ) , ( ˆ

1

i p i k p i p R

k i

          





 



) 1 , ( ) , ( ˆ

1

p p R p R

i i

 



 

• For i¥1 one has

) , ( ˆ ) , (    p R i p p R

i i



Interpolatory vector functions p y

• The polynomials with interpolating nodes as shown in figure

are of global order p=3 and have the form: are of global order p 3, and have the form:

1 2 1 1 ) , 2 ( ˆ ) , 2 ( ˆ ) , 2 (

3 2 1

    k j i p R p R p R

k j i

   2 1 , , 2 , 1 , ; , , 1 ,        p k j i p k j p i with  

Interpolation nodes for curl- or divergence-conforming bases on triangular elements Only nodes in basis subset Ω1

• r Λ1

for p = 3 are shown

• elements. Only nodes in basis subset Ω1

i j k or Λ1 i j k for p = 3 are shown.

82

References References

• R.D. Graglia, D.R. Wilton and A.F. Peterson, “Higher order

g g interpolatory vector bases for computational electromagnetics,” invited paper, special issue on “Advanced Numerical Techniques in Electromagnetics” IEEE Trans. Antennas Propagat., vol. 45, no. 3, pp. 329-342, Mar. 1997.

• R.D. Graglia, D.R. Wilton, A.F. Peterson, and I.-L. Gheorma, “Higher
• rder interpolatory vector bases on prism elements,” IEEE Trans. Antennas

Propagat vol 46 no 3 pp 442-450 Mar 1998 Propagat., vol. 46, no. 3, pp. 442-450, Mar. 1998.

• R. D. Graglia, and I.-L. Gheorma, “Higher order interpolatory vector

bases on pyramidal elements,” IEEE Trans. Antennas Propagat., vol. 47,

• no. 5, May 1999.

83

Volumetric Elements Volumetric Elements

84

Volumetric Elements

85

Volumetric Elements Volumetric Elements

86

Volumetric Elements

87

Example of results for surface elements:

3

INHOMOGENEOUSLY FILLED WAVEGUIDES



2.5 5

b h r 0

1.5 2 kz/k0 m=1 2 3 4 3 4 5

a=2b, h=0.1b a

0.5 1 6 7 8 1 1 2 6

r =10

1 2 3 4 5 6 7 8 9 10 k0*a 3 2

analytical + * FEM

88

+, FEM

• INHOM. FILLED WG - RELATIVE ERROR

10 10

• 1

OR

r 0 P 1

10

• 3

10

• 2

LATIVE ERRO

P=1

10

• 4

10 REL

P=2 P=3

10

2

10

3

10

4

10

• 5

MATRIX DIMENSIONS

89

Normal field component at the air-dielectric interface Normal field component at the air dielectric interface.

Fundamental mode

15

P=2 P=3

a b h r 0

Fundamental mode

10

ya/Eyd

P=1

a

a =2b, h =0.2 b 10 k 7

5

P 0

~ 1800 UNKNOWNS

Ey

P=1

r =10, k0a =7

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

P=0

x

P N 1641 1 1873 90 1 1873 2 1897 3 1681

Image waveguide

40

a

• - p=0, (37 incognite), mesh densa
• p=3, (325 incognite), mesh lasca

espansione modale [1]

35 40 z

b a` b`

a =1.3 mm, b=1.6 mm, a` =0.55 mm, b`=0.82 mm rx =170, ry = ry 85

3 4

30 ENZA GHz

modo 2

25 FREQUE

modo 1

2 4 6 8 20

1 [1] J I Askne E L Kolberg L Pettersson ``Propagation in a waveguide partially filled with

24 triangoli 19 nodi 291 triangoli 166 nodi 91

Kz mm-1

[1] J.I.Askne, E.L. Kolberg, L. Pettersson, Propagation in a waveguide partially filled with anisotropic dielectric material``, IEEE Trans. MTT, vol.30, n5, pp.795-799, Maggio 1982.

Singular vector bases Singular vector bases

92

References

[1] R.D. Graglia, G. Lombardi, “Singular Higher-Order Complete Vector Bases for Finite Methods,” IEEE Trans. Antennas Propagat., Vol. 52, No. 7, pp. 1672-1685, J l 2004 July 2004. [2] R.D. Graglia, G. Lombardi, “Hierarchical singular vector bases for the FEM solution of wedge problems ” invited paper for Proceedings of 2004 International solution of wedge problems, invited paper for Proceedings of 2004 International Symposium on Electromagnetic Theory, URSI-Commission B, Pisa, Italia, 23-27 May, 2004. [3] R.D. Graglia, G. Lombardi, “Singular Higher Order Divergence-Conforming Bases of Additive Kind and Moments Method Applications to 3D Sharp-Wedge Structures,” IEEE Trans. Antennas Propagat., vol. 56, no. 12, pp. 3768-3788, , p g , , , pp ,

• Dec. 2008. DOI 10.1109/TAP.2008.2007390

93

Circular vaned waveguide – FEM application =1/2 =0

94

element’s type and conformity

Circular vaned waveguide: eigenmodes =1/2

• Singular behavior
• Numerical precision
• Regular mode
• Singular mode

95

Modeling capability: very small thickness Double-vaned Circular homogeneous waveguide =2/3

96

Modeling capability: Multiple singular verteces and curvilinear singular elements curvilinear singular elements =2/3 =/2 = =1/2   2 /3 =3/4 =2/3

97

The square PEC-plate problem at normal incidence – MoM Application pp

The results at left (a c) were The results at left (a, c) were

• btained by using the zeroth-
• rder regular base (p = 0) on the

dense mesh A. The results at right (b, d) were

• btained by using the coarse

mesh B and the singular base of

98

• rder [p = 2, s = 0].

The square PEC-plate problem at normal incidence

99

(10 × 1) PEC-strip

Normal incidence Ex Singualar bases p = 2, s = 0 g p ,

100

The circular PEC-plate at normal incidence

101

The circular PEC-plate at normal incidence (d = /100)

102

Normal incidence on a (1 × 1) square PEC-plate with a hole of radius r = /10 centered at (x = −0.15, y = +0.15); the incident magnetic field is polarized in the y-direction.

103

Spherical PEC-shell of radius a = /(2 and aperture angle = 120° illuminated by a planewave propagating in the positiveˆz direction y p p p g g p

104