EP1290590A1

EP1290590A1 - Algorithme de simulation electromagnetique, notamment des performances d'une antenne

Info

Publication number: EP1290590A1
Application number: EP01940675A
Authority: EP
Inventors: H. Thales Intel. Property CHRISTIANSEN SNORRE; François Thales Intellectual Property BEREUX; Jean-Claude Thales Intellectual Property NEDELEC; Jean-Paul Thales Intelletual Property MARTINAUD
Original assignee: Thales SA
Current assignee: Thales SA
Priority date: 2000-06-09
Filing date: 2001-06-01
Publication date: 2003-03-12
Also published as: FR2810136B1; US20020198670A1; WO2001095159A1; FR2810136A1

Abstract

La présente invention se rapporte à un algorithme de simulation électromagnétique, qui permet de calculer l'onde électromagnétique diffractée par un conducteur en régime monofréquentiel. A cet effet, l'invention concerne un algorithme de simulation électromagnétique basé une résolution itérative d'un système d'équations intégrales comprenant un préconditionneur. Ce préconditionneur est notamment issu de l'adaptation des formules de Calderon aux équations intégrales de frontière de l'électromagnétisme, encore connues sous le nom de Electric Field Integral Equation (EFIE) dans la littérature anglo-saxonne. On utilise aussi une représentation origniale du résidu des calculs lors de chaque itération. Cette représentation, ainsi qu'une projection et une composition interviennent dans l'expression dudit préconditionneur. Elle s'applique notamment aux outils de simulation utilisés lors de la conception d'antennes de réception ou d'émission telle que les antennes de téléphone cellulaire, de radar anti-collision, de système de contre mesure électronique (CME), de radar de veille ou de poursuite, ou de satellite. L'invention s'applique aussi au calcul des surfaces équivalentes radar (SER) d'objets dont on connaît les propriétés géométriques.

Description

ALGORITHME DE SIMULATION ELECTROMAGNETIQUE, NOTAMMENT DES PERFOR MANCES D'UNE ANTENNE.

La présente invention se rapporte à un algorithme de simulation électromagnétique, notamment des performances d'une antenne, qui permet de calculer l'onde électomagnétique diffractée par un conducteur en régime monofréquentiel. Elle s'applique notamment aux outils de simulation utilisés lors de la conception d'antennes de réception ou d'émission telle que les antennes de téléphone cellulaire, de radar anti-collision, de système de contre mesure électronique (CME), de radar de veille ou de poursuite, ou de satellite. L'invention s'applique aussi au calcul des surfaces équivalentes radar (SER) d'objets dont on connaît les propriétés géométriques. On utilise les simulations d'antennes pour limiter le nombre de maquettes et de prototypes lors de la conception desdites antennes. Ces simulations permettent notamment de calculer le diagramme de rayonnement en champ lointain des antennes et d'adapter les antennes en émission ou en réception, en présence ou non d'une structure environnante. Elles utilisent comme donnée d'entrée un maillage de l'antenne dont on souhaite évaluer les performances, ainsi que les caractéristiques de l'excitation électromagnétique à laquelle elle est soumise. L'invention ne se limite pas aux simulations d'antennes. Elle s'applique aussi par exemple aux calculs de SER de cibles. On décrit maintenant, à titre d'illustration, l'application de l'invention dans les simulations d'antennes, en réception.

On distingue deux principales méthodes dans les simulations couramment employées. Une première méthode est basée sur le calcul par différences finies, encore connue sous le nom de méthode d'éléments finis de volume. Selon cette méthode, on utilise un maillage d'un volume entourant l'antenne. Un inconvénient de cette méthode est que le maillage est nécessairement borné, alors qu'on s'intéresse au diagramme de rayonnement à l'infini. On doit alors réaliser un compromis entre la dimension du volume maillé, c'est à dire la précision de calcul, et le temps de calcul. Pour pallier cet inconvénient, on utilise une seconde méthode basée sur des équations intégrales dans le domaine fréquentiel. Selon cette méthode, on utilise un maillage en surface de l'antenne uniquement. On calcule directement le diagramme de rayonnement à l'infini à partir de courants électriques et magnétiques sur la surface de l'antenne.

Certaines techniques connues utilisant les équations intégrales calculent les courants électriques et magnétiques (desquels on déduit le diagramme de champ rayonné à l'infini) grâce à une factorisation d'une matrice. Cette matrice est connue sous le nom de matrice d'interaction ou encore matrice d'impédance. Cette factorisation permet un calcul direct, c'est à dire non itératif des courants de surface. Un inconvénient de ces techniques est que le temps de calcul est long. Si on note N le nombre de points intervenant dans le maillage de l'antenne (encore appelé triangulation de surface), le temps de calcul selon ces techniques varie comme N³. Or le nombre de points N est lui-même lié à la longueur d'onde (et par conséquent la fréquence) de l'onde rayonnée par l'antenne. Supposons que l'on réalise une simulation à 10GHz, en utilisant N points dans le maillage de l'antenne, et que le temps de calcul est τ. Pour transposer cette simulation à 20GHz, il sera nécessaire d'utiliser un maillage comportant 4xN points, ce qui représentera un temps de calcul de l'ordre de 4³xτ. Un problème de temps de calcul se pose aussi lorsqu'on cherche à simuler des géométries complexes d'antennes, telles que les réseaux de petites tailles. Ceci rend ces techniques inutilisables notamment dans les outils de conception qui requièrent un temps de calcul réduit pour permettre aux concepteurs de réaliser plusieurs tests.

D'autres techniques connues utilisant les équations intégrales permettent de diminuer le temps de calcul grâce à une méthode de résolution itérative. Si on note IT le nombre d'itérations, le temps de calcul selon ces techniques varie comme ITxN². Un problème de ces techniques est que rien ne garantit la convergence des calculs. En d'autres termes, il existe des formes d'antennes pour lesquelles on ne peut pas calculer le diagramme de champ rayonné avec ces techniques.

Un but de l'invention est de pallier les inconvénients précités, et notamment de réduire les temps de calcul.

A cet effet, l'invention concerne un algorithme de simulation des performances d'une antenne basé sur une résolution itérative d'un système d'équations intégrales comprenant un préconditionneur. Ce préconditionneur est notamment issu de l'adaptation des formules de Calderon aux équations intégrales de frontière de l'électromagnétisme. En particulier dans le cas d'une antenne entièrement métallique, on propose un préconditionneur pour l'équation appelée « Electric Field Intégral Equation » (EFIE) dans la littérature anglo-saxonne. On utilise aussi une représentation originale du résidu des calculs lors de chaque itération. Cette représentation, ainsi qu'une projection et une composition interviennent dans l'expression dudit préconditionneur.

L'invention a les principaux avantages suivants : " elle converge rapidement ;

• elle permet de simuler des géométries arbitraires et des excitations quelconques ;

• son conditionnement est indépendant de la finesse du maillage ;

• elle s'adapte aux algorithmes basés sur le calcul d'une matrice d'impédance, en réutilisant ladite matrice d'impédance ;

- elle permet de traiter des antennes comportant, en plus du métal, des matériaux diélectriques.

D'autres caractéristiques et avantages de l'invention apparaîtront plus clairement dans la description qui va suivre et dans les figures annexées qui représentent :

- la figure 1 , un maillage d'une antenne ;

- la figure 2, une vue en coupe du maillage de la figure 1 ;

- la figure 3, un détail du maillage de la figure 1 sur lequel est représenté un champ de vecteurs ; - la figure 4, un diagramme fonctionnel d'un algorithme itératif ;

- la figure 5, une fonction de base constante par triangles sur un maillage ;

- la figure 6, une fonction de base affine par triangles et continue sur un maillage ; - les figures 7 et 8, deux illustrations des performances de l'algorithme selon l'invention comparé aux techniques connues.

On se réfère maintenant aux figures 1 et 2 qui représentent un exemple de forme d'antenne dont on cherche à déterminer le champ diffracté lorsque l'antenne est éclairée par une onde incidente. En d'autres termes, on cherche à simuler une antenne en réception. L'antenne prise dans cet exemple est une cavité sphérique de demi-angle d'ouverture π/4. Le rayon intérieur est 7/8, le rayon extérieur est 9/8 (unité de longueur arbitraire). La surface de l'antenne est notée r. Cette surface r est maillée par des triangles. La surface r est supposée être celle d'un conducteur parfait (l'antenne) Ω- plongé dans le vide Ω+.

Dans cet exemple, l'onde électromagnétique incidente qui illumine l'antenne est une onde plane monofréquentielle. Cette onde électromagnétique incidente, de nombre d'onde k^ιπc connu, est représentée par deux champs de vecteurs notés E^ιnc et H^ιnc correspondant respectivement aux champ électrique et au champ magnétique. Bien entendu, l'invention ne s'applique pas aux seulement aux ondes planes. Il est possible de substituer cette onde incidente par le champ émis par un dipôle rayonnant par exemple (mode émetteur de l'antenne).

On cherche à déterminer l'onde électromagnétique diffractée par l'antenne à l'infini, c'est à dire le diagramme de rayonnement en champ lointain. Cette onde électromagnétique diffractée est représentée par deux champs de vecteurs notés E^dlf et H^dlf correspondant respectivement aux champ électrique et au champ magnétique.

A partir du champ de courant circulant à la surface de ladite antenne, encore appelé courants électriques et magnétiques de surface, il est possible de calculer le champ électromagnétique rayonné en tout point de l'espace. L'expression asymptotique à l'infini du champ électromagnétique rayonné est le diagramme de rayonnement que l'on cherche à déterminer. Ce calcul bien connu en électromagnétisme est rappelé dans le document « Intégral Equation Methods In Scattering Theory » de D. Colton et R. Kress - John Wiley & Sons, New-York, 1983.

Le champ électromagnétique satisfait les équations de Maxwell dans le vide Ω+ qui s'écrivent :

rôt(Ë)= iωμH (1)

rôt(H)= -iωεË (2)

Les termes utilisés dans ces équations représentent : rot l'opérateur rotationnel ;

E = E^,πc + E^dlf le champ électrique total en notation complexe ; H = H^ιπc + H^dlf le champ magnétique total en notation complexe ; • ω la pulsation de l'onde électromagnétique ;

• μ la perméabilité magnétique ; ^» ε la permittivité électrique.

On rappelle que le nombre d'onde est relié à la pulsation 5 simplement par la relation suivante :

k = ω^με (3)

Le champ électrique de l'onde électromagnétique diffractée o E^dlf s'exprime en tout point B de Ω+ à partir des courants de surface notés u par les relations suivantes :

Ë^dif(B) = k |G_k(A_IB) û(A) dS_A +-lgrad_B M | JG_k(A,B) div_A(û) dS (4)

Aer A ASs:r

1 e ik AB

I5 G_k(A,B) = (5)

4π AB

Les termes des relations (4) et (5) représentent : div_A(ϋ) la divergence prise en un point A du champ de vecteur ϋ ; grad_B le gradient pris en un point B ; 0 » dS_A un élément différentiel de surface ; G_k la fonction de Green standard; AB la distance entre les points A et B ; k la norme du nombre d'onde défini par la relation (3).

L'homme de l'art pourra aller chercher d'autres éléments 5 techniques relatifs à ce calcul dans le document susmentionné dans la mesure où celui-ci fait partie intégrante de la description.

Le calcul des courants de surface, c'est à dire du champ de vecteurs û , est déterminé à partir de l'équation variationnelle suivante : 0

Vv m(ϋ,v) = l(v) (6) dans laquelle

m(ϋN) = k jjG_k(A,B) û(A)- v(B) dS_A dS_E

Asr B≤Γ

A jjG_k(A,B) div_A(ϋ) div_B(v) dS_A dS_E (7)

K Aer

Ber

et

l(v) = - JË^inc(A). v(A) dS_y (8)

Cette équation variationnelle (6) est connue sous le nom de d'équation intégrale de frontière de l'électromagnétisme, ou encore « Electric Field Intégral Equation » (EFIE) dans la littérature anglo-saxonne.

Afin de résoudre cette équation variationnelle (6) dans une simulation numérique d'antenne, on doit approcher la solution û dans un espace de dimension finie, dit espace de discrétisation. Cet espace contient des champs de vecteurs qui représentent des courants de surface. La dimension cet espace est le nombre de composantes nécessaires pour décrire entièrement ledit champ de vecteurs, en tout point de la surface r. Cette surface étant maillée, le nombre de composantes servant à décrire ϋ sera fonction notamment du nombre de points du maillage N, ainsi que de la nature des fonctions de base servant à décrire le champ de vecteurs (par exemple des fonctions linéaires ou de degré 2). Dans la suite de la description, on prendra à titre d'illustration l'espace de Raviart-Thomas de plus bas degré sur le maillage de r. Cet espace de discrétisation, qui représente des courants de surface, est noté Dh. Le symbole h représente la longueur caractéristique, ou encore précision, du maillage. En effet, le nombre de dimension de cet espace Dh dépend du nombre de point N du maillage, qui lui-même dépend de la précision h du maillage. On se réfère maintenant à la figure 3 pour décrire une base de l'espace Dh. Dans cette base, le champ de courants de surface ϋ est représenté par les coefficients d'un vecteur noté U.

La base de l'espace des courants de surface contient des éléments notés φ_i où i est un indice entier associé à une arête du maillage de la surface r. Ces éléments sont des champs de vecteurs définis sur le maillage de la surface r. Le champ de vecteurs φ_{i (} représenté par des flèches sur la figure 3, a un support borné à deux triangles Ti et T₂ du maillage. Ces triangles Ti et T₂ ont en commun l'arête d'indice i et de longueur £_s . Cette arrête est orientée par un vecteur de norme unitaire noté

Y_j . On note Pi le sommet du triangle Ti non contenu sur l'arête i ; on note P₂ le sommet du triangle T₂ non contenu sur l'arête i. On note Si et S₂ les surfaces des triangles Ti et T₂. Soient z₁ et z₂les vecteurs de norme unitaire, ayant une direction normale à la surface des triangles Ti et T₂, et orientés de l'intérieur Ω- vers l'extérieur Ω+. On défini les vecteurs x!, et

X₂ les vecteurs de normes unitaires tels que le triplet x^y^z et le triplet

(x₂;y_j;z₂ soient des trièdres directs. On définit le champ de vecteurs φ, pour tout point A appartenant à la surface T avec les relations suivantes :

- si A e ^, φ_i(A) = ±-^LPJA avec φ_i(A)- xⁱ ₁ > 0 ; (9)

- si A e T₂, c i(A) = avec φ.(A) - x₂ > 0 ; (10)

- sinon, φ, (A) = 5. (11)

Une telle base de l'espace des courants est connue sous le nom de base usuelle de l'espace Raviart-Thomas ou encore courants élémentaires de Rao-Wilton-Glisson. L'homme de l'art pourra aller chercher d'autres éléments techniques dans « Electromagnetic scattering by surfaces of arbitrary shape » de S.S.M. Rao, D.R. Wilton et A.W. Glisson - IEEE Trans. Ant. Prop. AP-30, pp. 409-418, 1982 - dans la mesure où ce document fait partie intégrante de la description. On note Uj les coefficients du vecteur U. Le champ de vecteurs u(A) utilisé dans la relation (4) se décompose sur la base usuelle de l'espace de Raviart-Thomas de la manière suivante :

û(A) = ∑U, φ,(A) (12)

On transpose maintenant en notations matricielles les relations (6), (7) et (8) en utilisant la base (cp_j) décrite ci-dessus. L'équation variationnelle (6) s'écrit sous la forme du système d'équations linéaires suivant :

M U = L (13)

Les termes de la relation (13) sont les suivants : ^» M est une matrice connue, dite matrice d'interaction ou encore matrice d'impédance ; " L est un vecteur connu, dont les coefficients représentent l'onde incidente, c'est à dire l'excitation électromagnétique ; • U est le vecteur que l'on cherche à déterminer, dont les coefficients représentent les courants de surface.

On définit les coefficients My de la matrice d'interaction M et les coefficients Lj du vecteur L représentant l'onde incidente par les relations suivantes :

M, = k JjG_k(A,B) (φ,(A). φ_J(B))dS_A dS_B

Aer Bεr

-Z JI^G _k(A.B) div_A(φ_i)diV_B(φ_j)dS_A dS_B (14)

K Aer Ber

L, = j^^»(A). φ,(A))dS_A (15)

Aer L'homme de l'art pourra aller chercher d'autres éléments techniques relatifs au calcul de ces coefficients dans « Approximation par éléments finis de surface de problèmes de diffraction des ondes électromagnétiques » de A. Bendali - Thèse de l'Université Paris VI, 1984 - dans la mesure où ce document fait partie intégrante de la description.

On se réfère maintenant à la figure 4 où est illustré un algorithme itératif de résolution du système d'équations linéaires (13) basé sur la technique de gradient conjugué. Il est à noter que la technique de préconditionnement que nous illustrons dans le cas de l'algorithme du gradient conjugué s'applique aussi bien d'autres algorithmes itératifs. On peut citer notamment les algorithmes du Résidu Minimum Généralisé, Gradient Bi-Conjugué, Résidu Quasi-Minimum et du Gradient Bi-Conjugué Stabilisé encore connus dans la littérature anglo-saxonne sous les noms de « Generalized Minimum Residual » (GMRES), « BiConjuguatè Gradient » (BiCG), « Quasi-Minimal Residual » (QMR), et « BiConjugate Gradient Stabilized » (Bi-CGSTAB). L'homme de l'art pourra aller chercher d'autres éléments techniques sur les méthodes itératives dans « Templates for the Solution of Linear Systems : Building Blocks for Itérative Methods » de R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donate, J. Dongarra, V. Eijkhout, R. Pozi, C. Romine et H. Van der Vorst - SIAM (1994), Philadelphia, PA - dans la mesure où ce document fait partie intégrante de la description.

La figure 4 illustre un algorithme 40 qui prend en entrée la matrice M et le vecteur L de l'équation (13) et donne en sortie le vecteur U. L'homme de l'art pourra aller chercher d'autres éléments techniques relatifs à la résolution de systèmes d'équations linéaires dans « Itérative Methods for Linear and Nonlinear Equations » de C.T. Kelley - SIAM Frontiers in Applied Mathematics, Philadelphia, 1995 - dans la mesure où ce document fait partie intégrante de la description.

Une première étape d'initialisation 41 permet d'initialiser quatre suites de vecteurs notés U[n], R[n], S[n] et P[n] où n est un entier. La première suite U[n] est une solution approchée qui converge vers la solution recherchée U. La seconde suite R[n], connue sous le nom de résidu, converge vers le vecteur nul. La troisième suite S[n], que l'on nomme résidu préconditionné, converge aussi vers le vecteur nul. La dernière suite P[n] est connue sous le nom de direction de recherche. Les premières valeurs de ces suites sont définies par les relations suivantes :

U[θ] = 0 (16)

R[0] = L - M U[0] (17)

S[0] = Z R[0] (18)

P[0] = S[0] (19)

La relation (16) est utilisée par défaut lorsqu'on ne connaît pas de solution approchée. Une variante de cette relation consiste à prendre pour U[0] le résultat d'un calcul de courant de surface réalisé pour une même antenne mais à une autre fréquence.

Le terme Z de la relation (18) est une matrice. Cette matrice est un préconditionneur pour la matrice M selon l'invention. On rappelle qu'un préconditionneur pour M est une matrice approchant l'inverse de M. Selon une variante, on peut remplacer le préconditionneur Z par l'identité. En d'autres termes, on peut initialiser S[0] avec R[0].

Une seconde étape d'itération 42 permet de calculer les valeurs des suites précitées à un rang n+1 à partir des termes de rang n. Cette étape d'itération utilise les relations suivantes :

U[n + l] = U[n] + α P[n] (20)

R[n + l] = R[n]- α M P[n] (21)

S[n + l] = Z R[n + l] (22)

r .₁ r i (Rfn + 1] , S[n + 11) _r -,

P n + 1] = S n + 1 ]+ X _r L. ^J/P n (23) ^{L J L J} (R[n] , S[n]) ^{L J}

avec (R n , S[nl α = ^L _Γ ^{J L} _r (24)

(M P[n], P[n)

où ( , ) représente le produit scalaire complexe en notation matricielle.

Une dernière étape 43 réalise un test de convergence. Ce test peut s'exprimer par exemple par l'inégalité suivante :

R[n] Il

^{L J} " < η (25)

R[0]

où η est un seuil prédéterminé. En d'autres termes, on compare le résidu préconditionné normalisé par rapport au un seuil prédéterminé η. Lorsque l'inégalité (25) est vérifiée, le calcul s'achève et on prend U = U[n] comme solution. Dans le cas contraire, on incrémente la valeur de n et on retourne à l'étape 42.

Bien entendu, il est possible d'utiliser un autre test de convergence pour arrêter les itérations. On peut par exemple remplacer l'inégalité (25) par l'inégalité suivante :

L'algorithme du gradient conjugué dans les techniques actuelles est utilisé sans préconditionneur, c'est à dire avec la matrice Z égale à l'identité. L'invention consiste à utiliser un préconditionneur basé sur une généralisation d'une formule de Calderon. L'homme de l'art pourra aller chercher des éléments techniques sur cette formule de Calderon dans « Mathematical Methods in Electromagnetism, Linear Theory and Applications » de M. Cessenat - World Scientific Publishing Co., page 89, 1996 - dans la mesure où ce document fait partie intégrante de la description. Cette formule s'écrit :

Dans cette relation (27), 94., J et ^ sont trois opérateurs sur les champs de vecteurs tangents, c'est à dire des applications qui à un champ de vecteurs tangents associent un autre champ de vecteurs tangents, et I est l'application identité. J est l'opérateur de produit scalaire avec la normale à la surface de l'antenne, et 94 est l'opérateur associé à la matrice d'interaction dont on cherche un préconditionneur. Les opérateurs 94, 3 et ^ sont définis formellement par les relations suivantes :

( M Û){B) = k |G_k(A,B) ϋ(A) dS

Aer

+ (A,B) div_A(ϋ) (28)

(7 Û)(B) = U(B) Λ Z(B) (29)

( û)(B) = jgrad_B(G_k(A,B)) Λ U(A)dS_A (30)

Aer

Dans les relations (28), (29) et (30), 94 ϋ , J ϋ, ^ϋ et ϋ sont des champs de vecteurs tangents. L'indice t de la relation (28) représente la composante tangentielle du vecteur entre crochets. B est un point de la surface r. z(B) est un vecteur unitaire, normal à la surface r en B orienté vers l'extérieur.

La Demanderesse a constaté que l'opérateur J <^J <B_^ étant compact, c'est à dire négligeable, l'opérateur 4 J 94 J est approximativement un inverse de l'opérateur 94. Ensuite, les préconditionneur étant définis à une constante multiplicative près, il est possible d'éliminer la constante 4. Sachant que J* = - J, le préconditionneur selon l'invention est défini à partir de * 4. J. On rappelle que J* est l'opérateur adjoint de , c'est à dire l'opérateur vérifiant la relation suivante :

Le préconditionneur selon invention est une traduction matricielle de l'opérateur J* 94 J. C'est dans cette traduction matricielle qu'apparaît l'avantage d'utiliser l'opérateur J* 94J plutôt que l'opérateur J 4J. En effet, la traduction matricielle Z de l'opérateur J* 94 J est une matrice symétrique, ce qui est essentiel pour les algorithmes itératifs, tel que le gradient conjugué.

On décrit maintenant un exemple de préconditionneur selon l'invention qui permet d'accélérer la convergence de l'algorithme et aussi de rendre cet algorithme plus stable (on converge toujours vers la solution quelles que soient les conditions initiales). Ce préconditionneur est adapté aux problèmes électromagnétiques et tire parti de la structure du problème à résoudre.

On définit d'abord des espaces et leurs bases associées qui serviront dans la suite de la description. On a déjà défini la base ( x).

L'espace engendré par cette base, Dh, est un espace de champs de vecteurs tangents. En d'autres termes, Dh est un ensemble de fonctions qui à tout point de r associent un vecteur tangent à la surface maillée r. On définit aussi deux espaces de fonctions Ch et Sh. Les fonctions de ces espaces associent une valeur scalaire à chaque point de r.

Ch est l'ensemble des fonctions constantes par triangles, c'est à dire dont la valeur est constante pour tout point appartenant à un triangle donné. La base de Ch que l'on utilise par la suite est l'ensemble des fonctions notées ψ_j , qui valent 1 sur le triangle d'indice i et 0 en dehors. On se réfère à la figure 5 sur laquelle est représentée une fonction de base ψ, .

On a représenté par soucis de clarté un maillage 50 plan et régulier. Les valeurs de la fonction ψ_j sont représentées sur un axe 51 , perpendiculaire au maillage 50. La fonction ψ_j représentée sur cette figure vaut 1 sur le triangle i et 0 en dehors.

Sh est l'ensemble des fonctions affines par triangles et continues. Ces fonctions ont un gradient constant sur tout triangle donné. La base de Sh que l'on utilise par la suite est l'ensemble des fonctions notées θ_j, qui valent 1 au nœud d'indice i et 0 sur les autres nœuds. On se réfère à la figure 6 sur laquelle est représentée une fonction de base θ_j . On a représenté par soucis de clarté un maillage 60 plan et régulier. Les valeurs de la fonction θ-, sont représentées sur un axe 61, perpendiculaire au maillage 60. La fonction θ_j représentée sur cette figure est affine par triangles, et vaut 1 au nœud i. Cette fonction a un support (des valeurs non nulles) borné aux triangles 62, 63, 64, 65, 66, 67 qui ont un sommet confondu avec le nœud i.

L'homme de l'art pourra aller chercher des éléments techniques complémentaires dans « Handbook of Numerical Analysis Vol. Il, Finite Eléments Methods (Part 1) » de P. G. Ciarlet - Ed. J.L. Lions, North-Holland,

1991 - dans la mesure où ce document fait partie intégrante de la description.

On définit maintenant des matrices qui serviront dans la suite de la description du préconditionneur. Ces matrices sont basées sur les bases définies ci-avant. Elles s'expriment par les relations suivantes :

M1„ = Jφ_i(A). φ_j(A) dS_A (32)

Aer

M2, = Jψ_j(A)div_A(φ_i)dS_A (33)

Aer

où z(A) est un vecteur de norme unitaire, normal à la surface r au point A, et orienté vers l'extérieur.

M4_{ij =}= J_Ψj(A) θ,(A) dS_A (35)

Aer

M5_{ij =} Jθ_j(A) θ_i(A) dS_A (36)

Aer

On décrit maintenant les étapes qui permettent de calculer le résidu préconditionné à partir du résidu dans les relations (18) et (22). Ces relations utilisent le préconditionneur Z que l'on décrit maintenant. Dans la description qui va suivre, on donne une décomposition de cette matrice Z en produit de matrices. Les matrices intervenant dans ce produit sont des matrices creuses ou des inverses de matrices creuses. Ainsi, selon une variante avantageuse de l'invention, on utilisera préférentiellement les matrices creuses plutôt que directement la matrice Z, ce qui permet notamment de réduire la mémoire utilisée et le temps de calcul. On note maintenant R le résidu, et S le résidu préconditionné. R correspond respectivement à R[0] et à R[n+1] dans les relations (18) et (22). S correspond respectivement à S[0] et à S[n+1] dans les relations (18) et (22).

Une première étape consiste à utiliser une représentation mixte. Cette étape sera mieux comprise à l'aide des notations vectorielles qui seront ensuite transposées en notations matricielles. Le résidu R correspond à une forme bilinéaire sur l'espace Dh défini ci-avant. On note p cette forme bilinéaire. On représente p par un champ de vecteurs r1 dans Dh de divergence nulle et une fonction q1 dans Ch d'intégrale nulle. Le champ r1 et la fonction q1 sont définis par la relation suivante pour tout champ de vecteurs V dans Dh :

Jrî(A) • v(A) dS_A + Jql(A) div _A (v) dS _A = p(v) (38)

Aer Aer

Le champ de vecteurs r1 ayant une divergence nulle, on peut écrire pour toute fonction f dans Ch et d'intégrale nulle :

Jf(A) div_A(ri)dS_A = 0 (39)

Aer

Dans la suite, pour ne pas compliquer inutilement la présentation, on ne rappelle pas que les champs scalaires considérés sont d'intégrale nulle. Les relations (38) et (39) se transposent de façon matricielle comme suit :

Dans la relation (40), on a utilisé une notation matricielle par blocs, dans laquelle : " *M2 est la matrice transposée de M2 ;

• 0 est un bloc nul (vecteur ou matrice) ; - R1 est la représentation matricielle de r ;

• Q1 est la représentation matricielle de q1.

On rappelle qu'il n'est pas nécessaire d'inverser la matrice définie par blocs pour déterminer R1 et Q1 à partir de R. En effet, cette matrice est une matrice creuse, c'est à dire qui contient beaucoup de termes nuls. La résolution d'un tel système est bien connue de l'homme du métier. Elle est rappelée dans « Handbook of Numerical Analysis Vol. Il, Mixed and hybrid methods » pp 523-640 de J.E. Roberts et J.-M. Thomas - Ed. J.L. Lions, North-Holland, 1991. Ce document fait partie intégrante de la description.

Une seconde étape consiste à projeter li et q1. Ces projections se traduisent en notations vectorielles par les relations suivantes :

r2 = _Dh(rl Λ zj (41)

q2 = ^_Sh(ql) (42)

ou :

^» _Dh est l'opérateur de projection dans Dh ;

" _Sh est l'opérateur de projection dans Sh ;

^» z est un vecteur unitaire, normal à r et orienté vers l'extérieur.

Ces projections définies par les relations (41 ) et (42) se traduisent de façon matricielle par : Cm o λ R2^ f M3 R[

(43)

0 M5 Q2 M4 Q1

La résolution de la relation (43) est similaire à la résolution de la relation (40) dans la mesure où les matrices M1 et M5 sont des matrices creuses. On peut noter que cette résolution est encore plus facile que celle de la relation (40) dans la mesure où les matrices M1 et M5 sont de plus symétriques, définies positives et bien conditionnées. On détermine ainsi R2 et Q2 à partir de R1 et Q1.

Une troisième étape consiste à combiner r2 et q2 en utilisant la relation suivante :

r3 = r2 - rot(q2) (44)

Cette troisième étape se traduit de façon matricielle par

R3 = R2 - M6 Q2 (45)

On appelle j l'application vectorielle qui résulte de la composition des trois étapes décrites ci-avant. Cette application est définie par :

j(r) = r3 (46)

Cette relation (46) se transpose en notations matricielles par :

J R = R3 (47)

où on note J la matrice correspondant à l'application vectorielle j, la matrice J étant définie par le produit suivant :

Il n'est bien entendu pas nécessaire de calculer directement J, car la décomposition (48) en produit de matrices creuses et d'inverses de matrices creuses est plus simple à utiliser. En d'autres termes, la matrice J est calculée de manière implicite lors de son utilisation.

On note z l'application vectorielle qui correspond à la matrice du préconditionneur Z. Le préconditionneur z est défini à partir de l'application vectorielle j par :

z = jO m o j (49)

où :

• ] ^' est l'adjoint de j ;

• o est l'opérateur de composition d'applications ; « m est l'application linéaire Dh → Dh^* associée à la forme bilinéaire m.

La relation (42) se traduit de façon matricielle par :

Z=^l J M J (50)

On se réfère maintenant aux figures 7 et 8 sur lesquelles sont illustrées les performances d'un algorithme avec un préconditionneur selon l'invention par rapport aux techniques connues n'utilisant pas de préconditionneur. Sur la figure 7, la courbe 70 représente la fonction n en l'absence de préconditionneur. La courbe 71 représente la fonction n -> log 10 sHP en présence du préconditionneur

|s[o]| décrit ci-dessus.

Sur la figure 8, la courbe 80 représente la fonction

en l'absence d'un préconditionneur. La courbe 81 r représente la fonction n H log 10 W en présence du préconditionneur

|R[o]|. décrit ci-dessus.

On constate que si on prend un critère de convergence classique η = 10^"4 dans la relation (25), 50 itérations suffisent avec le préconditionneur. Avec un algorithme non préconditionné, on n'atteint pas une précision suffisante en 200 itérations.

Ces simulations numériques on permis de montrer que l'utilisation d'un préconditionneur selon l'invention permet d'accélérer et de stabiliser les algorithmes itératifs. Bien entendu, l'invention ne se limite pas à l'exemple utilisé pour la décrire. Il est possible notamment d'utiliser d'autres fonctions de base ou un maillage différent de ceux pris en exemple.

L'invention se généralise à tout autre discrétisation de l'espace Dh. Si on prend pour Dh un espace différent de l'espace de Raviart-Thomas, on remplace les espaces Ch et Sh respectivement par :

Ch = {div(u) | u e Dh}+ 1 (51)

Sh = p I rôt(p) e Dh j (52)

En d'autres termes, • Ch est un espace d'éléments finis minimal tel que la divergence des éléments de Dh soit comprise dans Ch ; " Sh est un espace d'éléments finis maximal tel que le rotationnel des éléments de Sh soit compris dans Dh.

Une application principale de l'invention se trouve dans les outils de conception d'antenne, mais l'invention ne se limite pas à cette seule application. L'invention s'applique bien entendu aussi à tout outil de simulation basé sur le calcul du champ rayonné par un conducteur. On peut citer notamment le calcul de surfaces équivalentes radar (SER) d'objets dont on connaît les propriétés géométriques. Il est à noter aussi que la technique de préconditionnement décrite dans les cas d'une méthode itérative s'applique aussi à d'autres méthodes numériques rapides. Ces méthodes rapides sont basées sur une résolution itérative, mais seuls les termes utiles aux produits matrice-vecteur sont calculés. Ainsi, on calcule de l'ordre de N x log(N) éléments de la matrice d'impédance, au lieu de N² éléments selon les techniques itératives classiques. En d'autres termes, on calcule la matrice d'impédance de manière implicite. L'utilisation du préconditionneur selon l'invention dans ces méthodes rapides se fait sans difficulté. Ces méthodes sont intéressantes pour les objets de grandes taille, c'est à dire pour N grand. C'est le cas notamment pour les simulations d'antennes dites sur structure. On peut citer notamment les Méthodes Multipôles Multiniveaux et les Méthodes Intégrales Adaptatives, connues aussi dans la littérature anglo-saxonne sous les noms de « Fast Multilevel Multipole Methods » (FMM) et de « Adaptative Intégral Method » (AIM). L'homme de l'art pourra aller chercher des éléments techniques sur :

" les méthodes rapides en général dans « Fast Solution Methods In Electromagnetics » de W. C. Chew, J.-M. Jin, C.-C. Lu, E. Michielssen, J.M. Song - IEEE Trans. on Antennas and Propagation, 45(3):533-543, Mars 1997;

• les Méthodes Multipôles Multiniveaux dans « Multilevel Fast Multipole Algorithm For Electromagnetic Scattering By Large Complex Objects » de J. M. Song, C.-C. Lu, W.C. Chew, S.W. Lee - IEEE Trans on Antennas and Propagation, 45(10):1488-1493, Octobre 1997 ; ^» les Méthodes Intégrales Adaptatives dans « AIM : Adaptative Intégral Method For Solving Large Scale Electromagnetic Scattering And Radiation Problems » de E. Bleszynski, M. Bleszynski, T. Jaroszewicz - Radio Science, 31 (5):1225-1251, 1996 ; dans la mesure où ces documents font partie intégrante de la description. L'invention s'étend sans difficulté aux objets (antennes ou cibles) comportant des matériaux diélectriques. Dans ce cas, on cherche des courants électriques et magnétiques équivalents sur chaque interface. La matrice d'interaction entre ces courants comporte des blocs diagonaux du même type que la matrice M décrite ci-dessus. Un préconditionneur pour la matrice d'interaction est donc obtenu en considérant la matrice diagonale par blocs, dont les blocs sont du type de la matrice Z décrite ci-dessus.

Claims

REVENDICATIONS

1. Algorithme de simulation électromagnétique, pour déterminer l'onde électomagnétique diffractée par un corps en régime monofréquentiel, 5 à partir d'un maillage dudit corps et de l'excitation électromagnétique, caractérisé en ce qu'il comporte au moins :

(a) une détermination d'une matrice M, dite matrice d'interaction, dont les coefficients sont déterminés à partir du maillage dudit corps ;

(b) une détermination d'un préconditionneur Z de la matrice M, ce 1 o préconditionneur étant la traduction matricielle de l'opérateur J* 94J , où 94 est l'opérateur associé à M, J l'opérateur produit vectoriel avec la normale à la surface dudit corps, et J* l'opérateur adjoint de J; (b) une détermination des courants qui circulent à la surface dudit corps, par un algorithme itératif de type gradient conjugué utilisant ledit 15 préconditionneur Z, l'algorithme itératif permettant de résoudre une équation, dite d'équation intégrale de frontière de l'électromagnétisme, s'écrivant sous forme matricielle de la manière suivante :

M U = L

20. où U est un vecteur que l'on cherche à déterminer, dont les coefficients représentent les courants de surface, et L est un vecteur connu, dont les coefficients représentent l'excitation électromagnétique ; (c) une détermination de l'onde diffractée par ledit corps, à partir desdits 25 courants de surface.

2. Algorithme de simulation électromagnétique selon la revendication 1 , caractérisé en ce que les coefficients du vecteur U sont exprimés dans la base usuelle de l'espace de Raviart-Thomas.

30

3. Algorithme de simulation électromagnétique selon l'une des revendications précédentes, caractérisé en ce que le préconditionneur Z est déterminé de manière implicite.

4. Algorithme de simulation électromagnétique selon l'une des revendications précédentes, caractérisé en ce que le préconditionneur Z est défini par la relation suivante :

Z=* J M J

où J est une traduction matricielle de l'opérateur J, et ^lJ la matrice transposée de J.

5. Algorithme de simulation électromagnétique selon la revendication 4, caractérisé en ce que la matrice J est définie par la relation suivante :

ou

M1, = j _jtAHMdS,

M2_{ij =} J_Ψj(A)div_A(φ_j)dS_A

M3_{ij =} l(φ_i{A)Λ z(A)y φ_i(A) 4S_/

Aer

OU <f>_j est la base usuelle de l'espace de Raviart-Thomas ;

θi est la base des fonctions affines par triangles ;

5 ψi est la base des fonctions constantes par triangles ;

A un point appartenant à la surface T dudit corps ;

z(A) est un vecteur de norme unitaire, normal à la surface dudit o coprs au point A, et orienté vers l'extérieur.

6. Algorithme de simulation électromagnétique selon l'une des revendications précédentes, caractérisé en ce que l'algorithme itératif utilisé est un algorithme rapide, du type méthode multipôles multiniveaux. 5

7. Algorithme de simulation électromagnétique selon l'une des revendications 1 à 5, caractérisé en ce que l'algorithme itératif utilisé est un algorithme rapide, du type méthode méthode intégrale adaptative.

0 8. Algorithme de simulation électromagnétique selon l'une des revendications précédentes, caractérisé en ce que le corps est une antenne dont on cherche à déterminer une forme optimale, en utilisant l'algorithme de simulation dans un outil de conception d'antennes.

5 9. Algorithme de simulation électromagnétique selon l'une des revendications 1 à 7, caractérisé en ce que le corps est un objet de forme connue dont on cherche à déterminer la surface équivalente radar (SER).

0

5