FR 2734379 A1 19961122 FR 9511377 A 19950928 273437g Procédé et machine pour simuler le comportement acoustique d'un corps entouré par un milieu fluide La présente invention concerne la méthode des éléments finis (MEF) et des variantes et des extensions de celle-ci. La méthode des éléments finis est une méthode numérique pour modéliser le comportement de systèmes physiques. Elle est généralement mise en oeuvre à l'aide d'un ordinateur. Plus précisément, l'invention concerne un procédé et une machine utilisant la méthode des élénents finis et des méthodes connexes pour simuler le comportement acoustique d'un corps entouré par un milieu fluide. Des simulations de ce genre sont utiles, par exemple, pour prédire la distribution d'ondes sonores qui sont diffusées ou rayon- nées par une structure entourée par de l'eau ou de l'air. Dans le domaine de l'acoustique structurale traitée par le calcul, la question de la modélisation efficace du champ acoustique dans de grands domaines extérieurs est restée un problème difficile depuis plus d'un quart de siècle. De nombreuses techniques ont été développées. D'un point de vue historique, celles auxquel- les on a accordé le plus d'attention dans les publication techniques sont les suivantes: la méthode de l'équation intégrale aux limites ou, comme on l'appelle souvent, la méthode des éléments limites (MEL), qui est basée sur l'équation intégrale de Helmholtz de surface appliquée à la surface de la structure; et la méthode des éléments infinis (MEI), qui est basée sur l'équation différentielle de Helmholtz (l'équation d'onde réduite), appliquée à des secteurs semi-infinis du domaine (qui sont les éléments infinis) qui sont extérieurs à une limite artificielle entourant la structure, avec des éléments finis entre la structure et la limite. La méthode des éléments limites est le procédé préféré pour la plupart des chercheurs et des concepteurs de programmes. Son avantage principal consiste en ce qu'elle est basée sur une formulation mathématiquement rigoureuse, c'est-à-dire la représentation intégrale de Helmholtz, qui (i) satisfait la condition de rayonnement de Sommerfeld et (ii) représente la solution exacte dans toute l'étendue du domaine extérieur. De plus, le fait qu'elle soit capable de réduire d'une unité la dimension- nalité du problème a longtemps été considéré comme un avantage de calcul important. Cependant, nous avons trouvé que cela est en réalité désavantageux, du fait de la beaucoup plus grande largeur de bande des équations, à cause de la connectivité globale inhérente du procédé. Au contraire, jusqu'à présent, on n'a jamais envisagé de représenter avec des éléments finis la solution exacte, ou même une solution approchée, à l'inté- rieur de l'élément lui-même. Sur la base d'hypothèses physiques ou de théories asymptotiques appropriées, les éléments finis étaient seulement destinés à fournir une condition approchée de non-réflexion sur une limite artificielle entourant la structure, permettant ainsi d'obtenir une solution suffisamment exacte dans la région finie à l'intérieur de la limite. La solution à l'exté- rieur de la limite est ensuite habituellement obtenue à partir de l'intégrale de Helmholtz extérieure, en utili- sant les données de pression et de vitesse calculées sur la surface de la structure. Deux types d'éléments infinis ont été utilisés précédemment pour des applications acoustiques: l'élément à décroissance exponentielle et l'élément "projeté". L'élément à décroissance exponentielle procure une approximation de la décroissance spatiale de l'ampli- tude de pression acoustique par une exponentielle, p3ce e, dans laquelle Y est un nombre positif ajusté de façon empirique et r est une coordonnée quelcon- que qui s'étend jusqu'à l'infini. Du fait que cette fonc- tion de décroissance n'est pas en accord avec la théorie acoustique, la précision se dégrade rapidement lorsqu'on s'éloigne de la limite artificielle. L'élément projeté est basé sur la formeasympto- tique du champ pour r grand, en utilisant les termes initiaux des fonctions de Hankel sphériques d'ordre - i 1 __2 -ikr inférieur, c'est-à-dire pc(- + 2 +...)e.Du fait r r de la projection, on prévoyait la possibilité d'utiliser l'élément avec n'importe quelle forme et orientation, c'est-à-dire que celui-ci n'était pas nécessairement aligné avec des surfaces de coordonnées sphériques, permettant ainsi aux côtés "radiaux" de différents éléments de partir de différentes origines sphériques. Cependant, nous avons trouvé que pour que cette représentation converge, (i) les éléments doivent se trouver à l'extérieur de la plus petite sphère entourant la structure, et (ii) les côtés de tous les éléments dans une mailledoivent se conformer à des lignes radiales partant de l'origine de la coordonnée r unique dans le développement. Du fait que la région de fluide qui se trouve entre la sphère circonscrite et la structure doit être maillée avec des éléments acoustiques de taille finie, la maille acoustique totale deviendrait très grande et aurait donc une très mauvaise efficacité, au point de vue du calcul, pour des structures ayant un rapport de forme élevé. Par conséquent, l'élément présente un intérêt pratique seulement pour les structures "massives", c'est- à-dire celles que l'on peut circonscrire étroitement par une sphère. Il s'aqit d'une limitation importante pour de nombreuses applications. Un problème supplémentaire de l'élément projeté consiste en ce que la génération d'éléments est relative- ment coûteuse. En outre, la génération des matrices pour l'élément projeté exige une inversion d'une matrice mal conditionnée. Cela occasionne des problèmes de quadrature à des fréquences très élevées. Dans la demande de brevet des E.U.A. en instance n 08/330160, intitulée "A 3-D Acoustic Infinite Element Based On A Prolate Spheroidal Multipole Expansion", cédée à la demanderesse de la présente demande, nous avons décrit un nouvel élément infini basé sur un développement multipolaire en coordonnées sphéroidales allongées. Nous avons trouvé que l'utilisation de cet élément infini conduit à une efficacité considérablement améliorée dans la modélisation acoustique de structures allongées. Cepen- dant, on n'obtient généralement pas la même efficacité lorsqu'on utilise un élément infini allongé pour modéliser une structure qui a une forme générale plate, semblable à un disque. Comme l'élément infini sphéroîdal allongé mentionné ci-dessus, l'élément infini sphéroîdal aplati qui est décrit ici rompt avec la tradition dans la mesure o il est basé sur une formulation mathématiquement rigoureuse, c'est-à-dire un nouveau développement multi- polaire (une série exponentielle en 1/r dans laquelle r est dans ce cas un rayon sphéroidal aplati) qui (i) satis- fait la condition de rayonnement de Sommerfeld et (ii) représente la solution exacte dans l'ensemble de l'exté- rieur du domaine, jusqu'à un sphéroïde aplati circonscrit minimal. Lorsqu'on incorpore dans un élément infini les termes initiaux du développement, c'est-à-dire un polynôme en 1/r multiplié par un facteur oscillatoire, on peut augmenter arbitrairement la précision sur toute l'étendue de l'élément, en augmentant le nombre de termes dans le polynôme, ce qui est analogue à une extension en p dans des éléments finis basés sur des polynômes, de type classique. L'élément sphéroïdal aplati procure une amélio- ration vis-à-vis de l'élément projeté en éliminant toutes les limitations pratiques précitées. (i) Il élimine la projection tridimensionnelle (3-D) qui mélange ensemble la direction infinie (gênante) avec les directions finies (commodes), (ii) il généralise la sphère circonscrite à un sphéroïde aplati circonscrit, afin d'accepter des struc- tures qui sont plates (rapport de forme élevé, dans lequel deux dimensions sont très supérieures à la troisième) aussi bien que massives (faible rapport de forme, c'est-à- dire que les trois dimensions sont toutes similaires), et (iii) il tire parti de coordonnées séparables pour réduire les intégrales en trois dimensions ou 3-D à des produits d'intégrales en deux dimensions ou 2-D et une dimension ou 1-D que 1 'on peut calculer de façon beaucoup plus économique, ce qui élimine simultanément le mauvais conditionnement. Nous pensons que la technique de l'invention est probablement la technique la plus efficace qui soit disponible à l'heure actuelle pour de grands problèmes d'acoustique structurale traités par le calcul. De plus, la technique qui est utilisée pour cet élément acoustique peut probablement être étendue à d'autres types de propagation d'ondes, comme en électromagnétisme et en élastodynamique. Notre nouvel élément s'intègre aisément dans n'importe quel programme d'éléments finis structu- raux. Notre nouvel élément infini est particulièrement utile lorsqu'il est complété par l'élément infini allongé mentionné ci-dessus. Un homme de l'art disposant de ces deux éléments infinis dans une seule bibliothèque d'outils de modélisation par le calcul pourrait aisément modéliser un comportement acoustique avec une efficacité élevée dans une gamme étendue de configurations géométriques de problèmes. Des programmes d'éléments finis, en général, sont bien connus et il n'est pas nécessaire de les décrire ici en détail. Il existe de nombreux ouvrages traitant de ce sujet, parmi lesquels le suivant: D.S. Burnett, Finite Element Analysis, Addison-Wesley, Reading, Mass., 1987 (réimprimé en 1988), qui est incorporé par référence. D'autres caractéristiques et avantages de l'invention seront mieux compris à la lecture de la description qui va suivre de modes de réalisation, donnés à titre d'exemples non limitatifs. La suite de la descrip- tion se réfère aux dessins annexés dans lesquels: La figure 1A est un organigramme illustrant les étapes fondamentales dans l'utilisation et le fonctionne- ment d'un programme d'éléments finis. La figure lB est un schéma synoptique d'un exemple de système pour mettre en oeuvre la méthode des éléments finis. La figure 2 représente une structure plate (semblable à un disque), de forme par ailleurs arbitraire, ainsi que sa plus petite sphère circonscrite S, ayant un rayon r0. Les figures 3a-3d représentent les surfaces de coordonnées fermées de quatre systèmes de coordonnées orthogonales. Ces surfaces sont: (a) une sphère; (b) un sphéroïde allongé; (c) un sphéroïde aplati; et (d) un ellipsoïde. La figure 4 représente un ensemble d'ellipses et d'hyperboles co-focales. Lorsqu'on les fait tourner autour de l'axe z, ces courbes forment respectivement les surfaces de coordonnées sphéroidale aplatie et hyperbo- loidale à une seule nappe du système de coordonnées sphéroidales aplaties. La figure 5 illustre les coordonnées sphéroida- les aplaties, r, G et 0, et des paramètres géométriques associés. Sur la figure, le rayon aplati r est défini par 1(c1 + c2). fc1 2 La figure 6 représente une structure plate (semblable à un disque), de forme par ailleurs arbitraire, et un sphéroïde aplati circonscrit minimal S, avec un rayon aplati r0. La figure 7 représente un élément infini sphéroïdal aplati illustratif. La figure 8 représente une maille caractéris- tique utilisant des éléments infinis conformes à l'inven- tion. La figure 9 montre la base d'un élément infini illustratif, s'étendant sur le sphéroïde d'éléments infinis. Les lignes en pointillés sur la figure sont des lignes de coordonnées Q, 0. La figure 10 est une coupe en deux dimensions à travers une maille en 3 dimensions pour un problème d'acous- tique structurale classique. La figure 11 représente la discontinuité géomé- trique à une interface entre des éléments adjacents de taille infinie et de taille finie. La figure 12 illustre une convention de numéro- tation de noeuds locaux pour le cas particulier d'un quadrilatère de Lagrange quadratique (n = 9) dans les directions angulaires et d'un quadrupôle (m = 3) dans la direction radiale (N = m x n = 27). Exposé général de la méthode des éléments finis La méthode des éléments finis (MEF) est une méthode, mise en oeuvre à l'aide d'un ordinateur, pour prédire le comportement d'un système physique par l'obten- tion de solutions numériques d'équations mathématiques qui décrivent le système et ses conditions d'application de charges. On peut considérer que l'utilisation de la méthode des éléments finis comprend trois phases: le pré- traitement, la résolution et le post-traitement. On va maintenant envisager ces phases de façon plus détaillée en se référant à la figure 1A. Dans la phase de pré-traitement, le domaine physique du problème est divisé en une configuration de sous-domaines de forme géométrique simple, que l'on appelle des "éléments". La configuration résultante est ce que l'on appelle une "maille". De plus, des données du problème, telles que des propriétés physiques, des charges et des conditions aux limites sont spécifiées. Cette procédure est représentée par l'étape 10 de la figure. La phase de résolution comprend les étapes 20-50 sur la figure. Dans la première de ces quatre étapes, des valeurs numériques sont calculées pour les coefficients dans "l'équation matricielle d'élément" associée à chaque élément dans la maille. L'équation matricielle d'élément, dont l'établissement est indiqué dans la case 15 de la figure, est un ensemble de formules mathématiques pouvant faire l'objet d'un calcul numérique, qui sont établies de manière théorique et sont mises en oeuvre dans le program- me, ou code, d'ordinateur qui forme le programme de la méthode des éléments finis. Pendant l'utilisation du programme de la méthode des éléments finis, chaque élément dans la maille accède au code pour ces formules, et les valeurs numériques sont calculées pour les coefficients dans les formules en utilisant les données géométriques et physiques qui sont associées à chaque élément. Nous fournissons, pour la première fois, un ensemble de telles formules (c'est-à-dire l'équation matricielle d'élément) qui est spécifique d'un élément infini acoustique sphéroîdale aplati. La procédure utilisée dans l'établissement de l'équation matricielle d'élément qui est décrite en détail ci-dessous, met en oeuvre les idées générales suivantes. La variable de champ inconnue, pour laquelle l'analyse par la méthode des éléments finis recherche une solution (telle que la pression acoustique), est représentée de façon approximative à l'intérieur de chaque élément sous la forme d'une somme finie de fonctions connues, que l'on appelle des fonctions de "forme". On choisit habituelle- ment des polynômes pour ces fonctions de forme. Il y a dans cette représentation des paramètres inconnus, que l'on appelle des "degrés de liberté", qui deviennent les 273437g nouvelles inconnues pour lesquelles l'analyse par la méthode des éléments finis trouvera des valeurs. Les degrés de liberté sont souvent les valeurs que prend la variable de champ inconnue à des points spécifiques dans l'élément, que l'on appelle des "noeuds". Grâce à cette représentation, lorsque des valeurs pour les degrés de liberté seront ultérieurement calculées à l'étape 50, on connaîtra alors une solution en tout point, de manière continue, dans toute l'étendue de chaque élément. Ceci est possible du fait qu'à ce stade, on connaîtra à la fois les fonctions de forme et les paramètres, et que ceux-ci, considérés ensemble, définissent la solution complète. La représentation de la variable de champ inconnue, en termes de fonctions de forme, est ensuite introduite dans les équations physiques régissant le système (qui sont de façon caractéristique des équations différentielles ou intégrales) qui expriment les lois physiques auxquelles le système physique est soumis. Ces équations de calcul se réduisent à un système d'équations algébriques linéaires qui, lorsqu'elles sont écrites sous forme matricielle, constituent l'équation matricielle d'élément. Les expressions pour les coefficients dans cette équation matricielle sont manipulées par le calcul et par l'algèbre jusqu'à ce qu'elles soient sous une forme qui se prête à l'évaluation numérique par l'ordinateur. Ces coefficients numériques sont l'information de sortie de l'étape 20. A l'étape 30, les équations matricielles d'élément provenant de tous les éléments sont combinées ensemble (on dit qu'elles sont "assemblées") en une grande "équation matricielle de système". Les matrices associées à des éléments qui se touchent dans la maille se chevau- cheront partiellement, ce qui établira la continuité de la variable de champ d'élément en élément. Du fait que le chevauchement est partiel, la matrice de système a des 273437g dimensions qui augmentent de plus en plus au fur et à mesure que l'on assemble des équations d'élément, en donnant ainsi une seule très grande matrice de système. L'équation matricielle de système est ensuite modifiée (étape 40) pour tenir compte des conditions aux limites et des conditions d'application de charges. L'équation matricielle de système est ensuite résolue (étape 50), en utilisant des techniques classiques d'analyse numérique. Bien qu'il puisse y avoir des milliers, ou même des centaines de milliers, de degrés de liberté inconnus, la résolution de cette équation matri- cielle est généralement un problème qui ne présente pas de difficulté particulière. Ceci vient du fait que les éléments de la matrice ont tendance à prendre des valeurs différentes de zéro seulement dans une bande relativement étroite le long de la diagonale de la matrice. Une mesure bien connue de la largeur de cette bande (et donc du temps exigé pour la résolution) est la "demi-largeur de bande en valeur effLicace", B rms Dans la phase de post-traitement, la solution de l'équation matricielle de système est visualisée sous une forme significative appropriée (étape 60). Dans cette phase, une autre information utile peut être élaborée à partir de la solution et visualisée de façon similaire. On présente ci-dessous une description détaillée de l'établissement de l'équation matricielle d'élément pour l'élément infini acoustique sphéroidal aplati de l'inven- tion. Des codes ou programmes d'éléments finis acous- tiques structuraux ont un large domaine d'applications utiles. En simulant le champ acoustique complet dans une structure et autour de celle-ci, et son interaction avec la structure vibrante, la conception de la structure peut être modifiée plus rapidement et plus efficacement (en comparaison avec les procédures de prototypage répétées, ll utilisées actuellement, qui sont coûteuses), pour amélio- rer ou optimiser les performances acoustiques, par exemple pour réduire le niveau de pression sonore global ou pour modifier le diagramme de directivité. Des applications importantes comprennent la réduction du bruit affectant l'environnement qui provient de machines, d'automobiles, d'aéronefs, etc., et la réduction du bruit émis par des produits de consommation bruyants, tels que des outils à moteur, des appareils ménagers et des ventilateurs d'ordi- nateurs personnels. Des équipements terminaux pour des systèmes de télécommunications (haut-parleurs, postes téléphoniques à haut-parleur, cabines téléphoniques publiques, téléphones cellulaires, etc.) peuvent être mieux conçus pour la maîtrise de l'écho, la discrimination vis- à-vis du bruit de fond, l'annulation active du bruit et le renforcement de la parole. Des produits dans l'industrie des loisirs peuvent être conçus de façon plus optimale pour obtenir une excellente fidélité sonore; ceci s'applique à des systèmes stéréophoniques, des instruments musicaux et la partie audio de systèmes multimédia domes- tiques. Il existe également des applications militaires concernant le sonar, la réduction du bruit produit par des sous-marins et les contre-mesures acoustiques. Système pour mettre en oeuvre la méthode des éléments finis En se référant à la figure lB, nous allons maintenir décrire un système utile pour mettre en oeuvre la méthode des éléments finis en général, et le procédé de l'invention en particulier. Un générateur de maille 100 exécute la phase de pré-traitement 10 de la figure lA. Il est de pratique courante d'utiliser comme générateur de maille un ordina- teur universel programmé. Des descriptions du système physique à modéliser peuvent être introduites dans le générateur de maille directement à partir d'un terminal d'utilisateur 102, ou bien elles peuvent être introduites à partir d'un dispositif de stockage de données 105, tel qu'un disque magnétique ou une mémoire numérique. Le géné- rateur de maille appellera de façon caractéristique une bibliothèque enregistrée, 107, d'éléments (c'est-à-dire de configurations nodales disposées à l'intérieur de secteurs ayantune section transversale qui a la forme d'un triangle et/ou d'un quadrilatère). La maille et d'autres données de sortie du générateur de maille sont enregistrées dans un dispositif de mémoire 110, qui est de façon caractéristique un disque magnétique et qui peut également être facultativement une mémoire d'ordinateur. La phase de résolution (étapes 20-50 de la figure 1A) est également accomplie de façon caractéristi- que dans un ordinateur programmé, représenté par l'élément sur la figure. Comme indiqué, une étape intermédiaire de la phase de résolution consiste dans l'évaluation des coefficients des équations matricieles d'élément. Ces coefficients sont enregistrés dans une mémoire numérique, qui est représentée par l'élément 120 sur la figure. Une propriété générale de la phase de résolution consiste en ce que le temps de calcul augmente de manière exponentielle lorsque la fréquence acoustique augmente. Pour cette raison, il sera souvent avantageux d'utiliser pour l'ordinateur, c'est-à-dire l'élément 115, un proces- seur parallèle, tel qu'un processeur massivement parallèle ou un processeur parallèle pouvant être développé. Le Tableau 1 donne une indication de la capacité que nécessite l'élément 120. Pour la modélisation acous- tique, chaque élément infini exige deux matrices: une matrice de raideur acoustique et une matrice de masse acoustique. Chacune a les mêmes exigences de mémoire, qui sont indiquées dans le tableau conformément au type d'élément infini. Chaque article du tableau est évalué en multipliant 16 octets (pour un mot en double précision) par le carré du nombre de degrés de liberté. A l'étape 30 de la figure 1A, les coefficients de matrice d'élément sont assemblés en une équation matricielle de système. L'élément de mémoire pour enregis- trer les coefficients de matrice de système est représenté par l'élément 125 sur la figure lB. Le résultat de la phase de résolution est une description des valeurs que prend une variable de champ acoustique sur toute l'étendue de la maille. Cette infor- mation est soumise de façon caractéristique à un post- traitement dans un élément de post-traitement qui est désigné par la référence 130 sur la figure, et elle est émise vers un dispositif de mémoire 135 ou un dispositif de visualisation 140. TABLEAU 1 Type d'élément acoustique infini Exigences de mémoire (octets) Triangle quadratique x dipôle 2304 Empirique quadratique x dipôle 4096 Quadrilatère quadratique x dipôle 5184 Triangle quadratique x quadrupôle 5184 Triangle quadratique x octopôle 9216 Quadrilatère quadratique x 11664 quadrupôle Triangle quartique x dipôle 14400 Quadrilatère quadratique x octopôle 20736 Quadrilatère quartique x dipôle 40000 Triangle quartique x octopôle 57600 quadrilatère quartique x octopôle 160000 DESCRIPTION DETAILLEE A. Le développement classique à l'extérieur d'une sphère On considère une structure arbitraire immergée dans un fluide homogène infini et vibrant à une fréquence constante W (figure 2). On désigne par S une sphère de rayon minimal r0 qui circonscrit juste la structure, comme représenté sur la figure 2. Il est bien connu que la pression diffusée et/ou rayonnée p à l'extérieur de S peut être représentée par le développement multipolaire suivant en coordonnées sphériques r, 9, 0: 273437g e -i Fn (, 0; k) (1 r n= r dans lequel k est un nombre d'onde acoustique, les F sont n des fonctions différentiables sans points singuliers, et la série converge pour r > r0. Du fait que l'équation (1) est valide seulement à l'extérieur d'une sphère circonscrite, des éléments infinis basés sur l'équation (1) doivent se trouver seulement à l'extérieur d'une sphère circonscrite. La région entre la structure et la sphère doit donc être remplie avec des éléments acoustiques de taille finie. Ceci limite l'utilité de tels éléments infinis à des structures qui sont "massives", c'est-à-dire qui remplis- sent la majeure partie de l'espace à l'intérieur d'une sphère. Des structures longues et/ou plates nécessite- raient un nombre excessif d'éléments de taille finie (ou de degrés de liberté pour des éléments de type p), pour remplir l'espace. Par conséquent, un développement multi- polaire différent est nécessaire, et ce développement doit être approprié pour les structures non massives. Comme des développements multipolaires sont limités à un système de coordonnées spécifique (du fait qu'ils procurent une forme fonctionnelle explicite pour la coordonnée "radiale" qui s'étend jusqu'à l'infini), la première étape consiste à sélectionner un système de coordonées approprié. B. Systèmes de coordonnées pour différents rapports de forme Les surfaces de coordonnées fermées concernant respectivement quatre systèmes de coordonnées orthogonales sont représentées sur la figure 3. Ce sont toutes des ellipsoïdes, c'est-à-dire des surfaces quadriques pour lesquelles les trois sections transversales principales orthogonales sont des ellipses, ce terme englobant des cercles à titre 273437g forme dégénérée. Le cas le plus général, faisant interve- nir trois sections transversales elliptiques différentes, est le cas (d). On a donné aux autres des noms spéciaux: sphères, sphéroïdes allongés et sphéroïdes aplatis. Les sphéroïdes allongés et aplatis sont générés en faisant tourner des ellipses respectivement autour de leur grand axe et de leur petit axe. Des coordonnées sphéroïdales aplaties peuvent circonscrire étroitement une structure quelconque pour laquelle deux dimensions sont similaires et la troisième est comparable ou plus petite, c'est-à-dire des structures qui varient d'une forme "massive" à une forme "plate". Les applicationsauxquelles on s'intéresse ici concernent des structures plates (ou massives), pour lesquelles le système de coordonnées approprié est le système sphéroidal aplati. Par conséquent, on décrit ci- dessous l'élaboration d'un élément infirni sphéroidal aplati. C. Un nouveau développement à l'extérieur d'un sphéroïde aplati 1. Coordonnées sphéroidales aplaties La figure 4 montre une famille d'ellipses et d'hyperboles co-focales, qui sont orthogonales. Lorsqu'on fait tourner les courbes autour du petit axe des ellipses (l'axe z), elles deviennent des surfaces quadriques que l'on appelle respectivement des sphéroïdes aplatis et des hyperboloïdes à une seule nappe. Combinées avec des plans passant par l'axe z, ces surfaces orthogonales sont les surfaces de coordonnées qui constituent le système sphéroïdal aplati. La figure 5 définit les paramètres pertinents pour les coordonnées sphéroidales aplaties r, O, 0: c1 + c2 r = rayon aplati = 2 (2) G = angle polaire = complément du demi-angle du cône asymptotique de l'hyperboloide (0. G L -r) 0 = angle de circonférence, autour de l'axe z, identique à l'angle sphérique (0 2 2 2 f = rayon du cercle focal, x + y = f, z = 0 a = demi-grand axe = r b = demi-petit axe = r - f (r > f) f = excentricité du sphéroïde = - (4) r Des surfaces à r constant sont des sphéroïdes co-focaux, des surfaces à G constant sont des hyperboloi- des à une seule nappe co-focaux, et des surfaces à constant sont des plans contenant l'axe z. La transformation de coordonnées cartésiennes en coordonnées sphériques aplaties est la suivante: x = rsin Ocos., y= rsinOsino, z = \r2-f2 cosO (5) Il faut noter que des traitements classiques de fonctions d'onde sphéroidales donnent souvent la préfé- rence à un choix différent pour les coordonnées radiales aplaties, c'està-dire r = \T - f2. On pourra voir par exemple l'ouvrage suivant: C. Flammer, "Spheroidal Wave Functions", Stanford University Press, Stanford, Calif., 1957, pages 6-10. Le choix représenté par l'équation (5) est essentiel dans l'invention, du fait que la coordonnée radiale aplatie classique r conduira à un developpement multipolaire divergent (c'est-à-dire un developpement de la forme représentée par l'équation (6) ci-dessous). 2. Un developpement multipolaire sphéroidal aplati On considère la même structure que sur la figure 2, mais on désigne maintenant par S un sphéroïde aplati ayant un rayon aplati minimal r0 qui circonscrit juste la structure (figure 6). On définit ainsi une limite inférieure du fluide. La pression diffusée et/ou rayonnée p à l'extérieur de S peut être représentée par le dévelop- pement multipolaire suivant en coordonnées sphéroidales aplaties r, G, 0: e -i Gn(eP;k) = rn=O r (6) avec des propriétés similaires à celles de l'équation (1). La série dans l'équation (6) converge absolument et uni- formément en r, G et 0, dans toute région satisfaisant la condition r > r0 + E >r0. Il y a une formule de récursion à quatre termes qui définit les fonctions Gn, avec n 23, exprimées en fonction de GO, G1 et G2. L'équation (6) est une généralisation de l'équa- tion (1), englobant cette dernière à titre de cas spécial. Alors que l'équation (6) convient pour une gamme étendue de formes, variant depuis des formes massives jusqu'à des formes arbitrairement plates, l'équation (1) convient seulement pour des formes massives. D'un point de vue géométrique, des sphéroïdes aDlatis ont une gamme étendue de formes, mais des sphêre3 n'ont qu'une seule forme. On peut voir que l'équation (6) est appropriée pour des structures plates en considérant la différence de réponse radiale le long du grand axe et du petit axe des coordonnées sphéroldales aplaties. Dans le plan du grand axe (z = 0), la coordonnée radiale aplatie, r, est la distance à partir de l'origine, a; par conséquent, r est identique au rayon sphérique. Dans n'importe quelle autre direction, r n'est pas la distance à partir de l'origine. En particulier, le long du petit axe (axe z) b est la distance à partir de l'origine. Du fait que r = b2 + f, alors près de la structure, o b/f c.l, on a dr/db - b/f. Par conséquent, r change plus lentement que b. On peut apprécier ceci en réexaminant la figure 4; du fait que des sphéroïdes sont des surfaces à r constant, une décrois- sance en 1/rn est étalée sur de plus grandes distances radiales à G = O qu'à G = Tr/2. En résumé, la vitesse spatiale de décroissance dans la direction radiale dans l'équation (6) est semblable à l'étalement sphérique dans le plan du grand axe, mais près du côté d'une structure plate, la décroissance est plus lente, comme on pourrait le prévoir. Dans le champ lointain, a, b et r s'approchent de l'égalité, ce qui fait que la décroissance s'approche de l'étalement sphérique dans toutes les directions. D. Configuration géométrique d'un élément acoustique infini sphéroidal aplati multipôle variable 3-D L'élément infini est représenté sur la figure 7, une coupe transversale 2-D d'une maille 3-D caractéristique à l'extérieur d'un corps 3-D est représentée sur la figure 8. Une face de l'élément, appelée la "base", doit être attachée à une surface sphéroidale aplatie de rayon aplati r1 entourant la structure, que l'on appelle le "sphéroïde d'iélérment infini" ou "sphéroïde de limitation", et la base doit se confor- rmer à la forme de la surface sphéroidale aplatie. La base peut être un quadrilatère curviligne (représenté ici) ou un triangle; il n'est pas nécessaire qu'elle se conforme aux lignes de coordonnées G, X (figure 9). Ceci permet de construire sur le sphéroïde une frame 2-D colaunent générale, comprenant des quadrilatères curvilignes et des triangles de n'importe quelle forme et orientation. On peut utiliser sur la base n'importe quelle représentation d'élément fini, par exemple des polynômes de Lagrange, empiriques ou hiérarchiques, de n'importe quel degré. (La configuration nodale quadra- tique de Lagrange qui est représenté ici, et dans des figures ultérieures, est utilisée dans un but illustratif). Suivant un mode de réalisation, le sphéroïde de limitation peut coïncider avec la limite intérieure du fluide. Le sphéroïde d'élément infini détermine sans ambiguïté un rayon focal, f, qui définit à son tour un système de sphéroïdes et d'hyperboloides cofocaux (voir la figure 4) qui sont utilisés dans la construction du reste de l'élément. Les faces latérales de l'élément sont les lieux d'hyperboles co-focales émanant des côtés de la al + base. Un multipôle d'ordre m (p = +... + e) r r exige m couches de noeuds qui se trouvent sur des sphéroi- des co-focaux de rayons allongés r1, r2,..., rm. La configuration nodale est identique sur chaque sphéroïde. La valeur m = 2 correspond à un dipôle, m = 3 correspond à un quadrupôle (comme représenté sur la figure 7), m = 4 correspond à un octopôle, etc. Un élément "multipôle variable" contient un nombre variable de couches de noeuds, permettant à l'analyste d'affiner la maille en direction radiale dans la région infinie, en augmentant l'ordre du multipôle, m, de façon analogue à une extension en p pour des éléments de taille finie. Enfin, il y a une face sphéroidale extérieure, S, à un rayon aplati r, qui s'éloigne à l'infini dans le développement théorique. Par conséquent, dans l'établis- sement de l'équation matricielle d'élément, l'élément commence avec une taille finie r] r t À, ce qui fait que la condition de rayonnement de Sommerfeld peut être appli- quée à la face extérieure, et ensuite la limite est obtenue pour 1 -- '. Des coordonnées classiques (illustrées pour un quadrilatère sur la figure 9) sont définies sur la section transversale de l'élément, c'est- à-dire sur des sphéroïdes, par la projection: n 8(trn) = E v,) O(,n) v(,), v=1 v1 v t(7) dans laquelle n est le nombre de noeuds dans le quadri- * latère ou le triangle, Gu '., sont les coordonnées angulaires sphéroidales aplaties du 9-ième noeud, et Xva(,) sont des polynômes d'interpolation. (Selon une %arian--:, on pourrait utiliser des fonctions de raccorde- ment dans le cas o les éléments sont hiérarchiques.) Du fait que la base, les couches nodales inter- médiaires et la face extérieure se conforment à des sphéroïdes co-focaux et les faces latérales sont les lieux d'hyperboles co-focales, l'élément est un cylindre droit dans l'espace r, G, 0 (ou l'espace rt,) et les inté- grales 3-D se séparent en intégrales angulaires 2-D bien conditionnées et en intégrales infinies 1-D. Un avantage supplémentaire consiste en ce que la projection ci-dessus est limitée seulement aux deux dimensions angulaires finies, c'est-à-dire les coordonnées G, 0, sur des surfaces sphéroldales. La coordonnée radiale utilise une projection séparée. Ceci évite le mauvais conditionnement numérique associé à des projections 3-D qui mélangent la dimension infinie avec les deux dimensions finies. L'utilisation de ces éléments infinisest illus- trée sur la figure 10, qui montre une maille caractéristi- que pour un problème acoustique structural. Cette figure est une coupe transversale d'une maille 3-D de la structure et de la région de fluide infinie qui l'entoure. La structure présente une symétrie axiale et elle consiste en deux plaques parallèles minces jointes à leurs bords par une coquille annulaire mince avec une section transversale conique arrondie. Certaines structures internes se trouvent entre les plaques. Il y a deux couches d'éléments acoustiques de taille finie entre la structure et le sphéroïde d'élément infini, et ensuite une couche d'éléments infinis, repré- sentés avec des lignes en pointillés, à l'extérieur du sphéroïde. Les lignes en pointillés sont toutes des hyper- boles. La première couche d'éléments de taille finie est générée en projetant de la même distance vers l'extérieur, le long de normales à la coquille, tous les noeuds qui se trouvent sur la surface extérieure de la coquille, ce qui crée une surface parallèle à la coquille. La seconde couche projette vers l'extérieur, en direction du sphéroïde, le long de normales, les noeuds qui se trouvent sur la surface parallèle. Ceci est une procédure simple pour la génération de maille automatique pour des surfaces convexes; des surfaces concaves exigeraient un générateur de maille 3-D plus perfectionné. Il peut y avoir ou non une discontinuité (en coordonnées, ' et donc dans la variable dépendante, c'est-à-dire la pression) entre les éléments infinis et les éléments de taille finie adjacents, le long du sphéroïde d'élément infini (figure 11), en fonction de la manière selon laquelle la configuration géométrique pour les éléments de taille finie est générée. Si ces éléments de taille finie utilisent une projection polynomiale classique basée sur les coordonnées cartésiennes globales des noeuds, alors les éléments ne se conformeront pas exactement à la forme du sphéroïde. Cependant, la projec- tion de l'équation (7) pour l'élément infini définit un jeu de coordonnées, 1 qui s'étendent exactement sur le sphéroïde. Il y a donc une mince "tranche" d'espace, comprenant des vides et des chevauchements, entre les éléments fini et infini, et les coordonnées, M dans les deux éléments ne coïncident pas exactement sur leur inter- face. Ceci ne constitue pas un problème, du fait que l'erreur due à cette discontinuité géométrique converge vers zéro, au fur et à mesure que la maille est affinée en h ou p, à la même vitesse que l'erreur (bien connue) à la limite d'un domaine courbe lorsqu'on utilise des éléments classiques employant une projection polynomiale. En fait, des calculs numériques simples montrent que la séparation spatiale maximale entre des valeurs, m identiques sur les faces d'éléments fini et infini adjacents est de façon caractéristique inférieure de plusieurs ordres de grandeur à des distances internodales. (L'utilisation de fonctions de raccordement pour la configuration géométrique des éléments de taille finie éliminerait évidemment la discon- tinuité. E. Equations de physique régissant le système Le comportement harmonique en fonction du temps (ei) est régi par l'équation 3-D de Helmholtz, V2p + k2p = 0 (8) dans laquelle k est le nomore d'onde (= t'./c), c est la vitesse du son (= B/), B est le module de compressi- bilité, e est la densité et p est l'amplitude complexe de la pression diffusée et/ou rayonnée: piff totale pinC s'il y a diffusion ray totale P = p = potale s'il y a rayonnement diff+ray =totale inc (9) p - p - p s'il y ales deux Pour garantir l'unicité de la solution, la pression doit satisfaire la condition de rayonnement de Sommerfeld à la "limite" extérieure, à l'infini: a + ikp = o 1, r àT _r r(10) dans laquaeL '- o minuscule, qui se lit "petit o", signifie "d'ordre inférieur à". Le r dans l'équation (10) est un r sphérique. Cependant, des coordonnées sphéroidales aplaties s'approchent de coordonnées sphériques lorsque r ->oo, ce qui fait que l'on peut utiliser l'équation (10) avec des coordonnées sphéroidales aplaties à la limite, lorsque r --a,. Lorsqu'on applique l'équation (10) à la face extérieure de l'élément S, au rayon aplati r, lorsqu'il tend vers l'infini, l'équation (10) peut être satisfaite par un choix approprié de fonctions de forme. F. Représentation de la pression par éléments finis 1. Expression générale; numérotation des degrés de liberté La pression diffusée et/ou rayonnée est repré- sentée de la façon suivante: N p(Err) = ( 1f(,,r)Pj() J=I avec fVj(,N,r) = Wr(),x)v(r) v = 1,2.. .,n; = 1,2..m n x m = N (12) Les wa(,7) sont ici des fonctions de forme "angulaires"'qui interpolent p sur des surfaces sphéroîdales co-focales jusqu'à la surface d'élément infini, et les 8V(r) ont des fonctions de forme "radiales"qui interpolent p le long de rayons hyperboliques. La continuité C inter-élément est établie en imposant la propriété d'interpolation: V(vv = ô, (13) (r') = 84, (14) La convention de numérotation des noeuds locaux qui est actuellement oréférée consiste à commencer avec les noeuds se trouvant à la base de l'élément et à progresser ensuite radialement vers l'extérieur, une couche de noeuds à la fois. Ceci est résumé dans le Tableau I et est illustré sur la figure 12. 2. Fonctions de forme angulaires Les fonctions V(,) s,,it des polynômes 2-D (de type empirique, Lagrange ou hiérarchique) classiques. Par exemple, pour les éléments consistant en quadrilatères de Lagrange quadratiques qui sont repré- sentés sur les figures 7, 9 et 12, les fonctions de forme angulaires sont: Wv(,n) = ta()t"'(n) v = 1,2,...,9 (cc =1,2,3; c' =1,2,3) (15) avec 1 1 t(u) = u(u-I); z2(u) = I -; t3(u) = -u(u + 1) (15A) 2 2 1A Si les fonctions YV(%,n) sont également utilisées pour les fonctions de projection X%(t,) dans l'équation (7), alors l'élément est "isoparamétrique" dans les directions angulaires. 3. Fonctions de forme radiales Les fonctions Vf(r) utilisent une forme tronquée de la partie radiale du développement multipolaire dans l'équation (6), c'est- à-dire un développement multipolaire du m-ième ordre: )= e É 'l(kr) ' g = 1,2,...,.. m (m 2 2) (16) j.L1 (kr) ikr Le facteur de phase e P ne doit pas obligatoirement être inclus; s'il est omis, il apparaîtra simplement dans chaque h, dans l'équation (20) ci-dessous. Les facteurs k dans les dénominateurs ne sont également pas néces- saires; ils sont inclus seulement pour faire en sorte que le terme h, soit sans dimension. J4' Les coefficients h.., sont déterminés par l'exigence de continuité C inter-élément. L'application de l'équation (14) à l'équation (16) donne m systèmes de m équations algébriques linéaires: [HI][S] = [I] (17) avec Hop = h, (18) Sg, = (kr,)- 1, (19) et [I] est la matrice identité. On a donc, [H] = [S]-, (20) Cette procédure définit m couches de noeuds (avec n noeuds sur chaque couche) s'étendant sur des sphéroïdes de rayons aplatis r1, r2,..., rm (voir la figure 12). A titre d'illustration, on considère un élément dipole (m = 2), "(r) = e(r-") [hg, + j A = 1,2 (21) kr (kr)2 L'inversion d'un matrice [Si de 2 x 2 donne: [H] = [il h2j 1 r5krf k2rir2j (22) 21 h22 r2- rl -kr2 -k2rlrA La procédure des équations (16)-(20) est celle qui a été utilisée jusqu'à présent. Cependant, une formu- lation hiérarchique aurait les avantages habituels de facilité d'affinement de la maille (par extension en p), d'amélioration du conditionnement numérique et d'élimina- tion de tous les noeuds extérieurs au sphéroïde d'élément infini. Le dernier avantage est particulièrement important du fait qu'il éliminerait la nécessité de placer des noeuds sur des trajectoires hyperboliques, ce qui n'est pas une technique classique dans des générateurs de maille d'éléments finis. Pour effectuer la conversion en une formulation hiérarchique, les directions angulaires utili- seraient les fonctions de forme hiérarchiques 2-D classi- ques pour des quadrilatères et des triangles. La direction radiale utiliserait la projection s = 1 - 2r1/r, qui projette de façon linéaire l'intervalle 1/r e [1/r1,0) sur l'intervalle E1-1,1), et elle utiliserait ensuite les fonctions de forme hiérarchiques 1-D classiques dans, en excluant la fonction linéaire qui est égale à l'unité à l'infini, c'est-à-dire à 3 = 1. 4. Matrices de raideur, de masse et d'amortissement du rayonnement, et équation matricielle d'élément On peut trouver des détails mathématiques supplémentaires concernant l'établissement de ces expres- sions dans l'Annexe ci-après, ainsi que dans l'article de D.S. Burnett, "A 3-D Acoustic Infinite Element Based on a Prolate Spheroidal Multipole Expansion", J. Acoust. Soc. Am. 96(4), octobre 1994. 1. Expressions formelles On établit les matrices pour l'équation matri- cielle d'élément en partant d'un élément de taille finie, c'est-à-dire avec la face extérieure sur un sphéroïde de rayon aplati f (voir la figure 7), et en prenant ensuite la limite lorsque -,o. Ainsi, en appliquant à l'équa- tion (8), sur un seul élément, la méthode du résidu pondéré de Galerkin, on obtient: lim f f f (BV2p+c2pp)yidV=O, i=1,2,...,N, (23) 2 2 en utilisant k2 = /B et les fonctions de test 't i. Bien qu'il soit utile d'utiliser, pour des fonctions de test, les fonctions de base Yj qui sont utilisées dans l'équation (7) et sont représentées dans l'équation (23), ce choix n'est pas le seul. Ainsi, d'autres fonctions de test peuvent être utilisées, au moins dans certains cas, pour obtenir des avantages secondaires. Un exemple consis- te à utiliser N i' c'est-à-dire les conjugués complexes des ' i. Ce choix réduirait légèrement le temps de calcul des intégrales R. dans l'équation (112) ci-après. La première intégrale est convertie en une intégrale de surface et en une autre intégrale de volume, en utilisant l'identité ViV2p=V'(VIVp)- VVi.Vp, et le theorème de la divergence. En substituant l'équation (11) dans les intégrales de volume, on obtient l'équation matricielle d'élément suivante: ([K] - o2[M])(p) = {F), (24) dans laquelle la matrice de raideur [K], la matrice de masse |M] et le vecteur de gradient de pression {Fi sont respectivement: Kij = lim | f | BV17i. V'jdV, Mij = lim | | f pViyjdV, Fi = lim || B14in dS (25) --±}an L'intégrale de surface pour AF) est étendue sur la totalité de la limite de l'élément, S(e) Elle se décompose en deux intégrales: l'une sur la face exté- (e) () _ - e rieure S(e) et l'autre sur les faces restantes S(e) - S(e) Fi = liJ f,Bain dS + lim 2 B)i PdS. On considère la première intégrale dans l'équa- tion (26). Lorsque r -- o, les coordonnées sphéroidales aplaties tendent vers des coordonnées sphériques, ce qui fait que Zp/ n --> ?p/ r et dS -r sin OdGdO, en désignant par r, O, 0 des coordonnées sphériques. Pour évaluer p/ ar lorsque -, Do, on substitue les équations (11), (12) et (16) dans l'équation (10), ce qui donne une forme plus forte de la condition de Sommerfeld, + ikp = 0 r- (27) dans laquelle le O majuscule, qui se lit "grand O", signifie "au moins d'ordre équivalent à". On substitue àp/ ôr dans l'équation (27) à ôp/ n dans la première intégrale. Le terme O(1/r) n'apporte aucune contribution à l'intégrale, du fait que \fi est O(1/r) et dSocr lorsque --. Dans le terme ikp restant, on substitue l'équation (11) pour le terme p. Dans la seconde intégrale dans l'équation (26), on substitue le reste de la quantité de mouvement linéaire à -p/ n, c'est-à-dire -p/ Dn = c2 un, en désignant par u l'amplitude de la composante normale du déplacement n de la particule. En effectuant ces substitutions, l'équation (26) devient: (F} = - iw[C] (P} + (D}, (28) avec Cij = pc lim | I qivjdS Di= (opc)2 lim f f viundS. (29) -+ S{, _j ') La matrice [CI est la matrice d'amortissement de rayonnement, représentant la perte de rayonnement jusqu'à l'infini. Le vecteur {Dû permet de spécifier u sur les faces latérales ou inférieure de l'élément. En substituant l'équation (28) dans l'équation (24), on obtient, pour l'équation matricielle d'élément, ([K] + ico[C] - c2[MI){P} = (D), (30) dans laquelle, en résumant, les expressions formelles pour les matrices de raideur, de masse et d'amortissement de rayonnement sont respectivement, K;-=iimJ'f Vxî.V) v K ij= lim f f f pBVi. VjadV, A V'(,) C1 PC1ImI f (31) Cij = pc lim I I i ij dS.(31) Les D. sont égales à zéro dans virtuellement toutes les i applications pratiques. 2. Transformation d'intégrales; expression finales Les opérations mathématiques restantes transfor- ment les intégrales de l'équation (31) en expressions que l'on peut évaluer de façon numérique. (Les intégrales Di sont négligées du fait qu'elles sont égales à zéro dans virtuellement toutes les applications pratiques.) On peut trouver des détails mathématiques dans l'Annexe; on présente ci-après une brève description des principales étapes. * Transformation des intégrales de l'équation (31) en coordonnées sphéroldales aplaties r, G, g. Les éléments différentiels de volume et de surface sont dV = J drdGd0 et dS = J dGdd, en désignant V s par J le jacobien de volume et par J le jacobien de v s surface. Jv et J,S ainsi que l'opérateur de gradient, V, V S dans K.. dans l'équation (31), sont exprimés en coordon- nées sphéroidales aplaties. Ceci sépare les intégrales 3-D pour K1 et M. en produits d'intégrales "angulaires" 2-D sur G, 0, et d'intégrales "radiales" 1-D sur r, et sépare l'intégrale 2-D pour Cij en un produit d'une intégrale angulaire 2-D et d'une fonction de r. * Etablissement d'expressions finales pour les intégrales angulaires. On transforme les coordonnées G, 0 en coordon- nées, s locales, en utilisant la projection de coordon- nées de l'équation (7) [voir la figure 9]. Les intégrales bien définies résultantes peuvent être intégrées de façon numérique en utilisant les règles de Gauss de la manière classique pour les éléments finis. * Etablissement d'expressions finales pour les intégrales radiales pour Kij et Mij, et pour la fonction radiale pour C... 1] On substitue l'équation (16) dans chacune des intégrales radiales et dans la fonction radiale, et on effectue diverses opérations algébriques. Certaines des intégrales deviennent des transformées de Fourier en sinus ou en cosinus, bien définies, que l'on peut évaluer par des algorithmes standards disponibles dans de nombreux progiciels mathématiques. Les autres intégrales, ainsi que la fonction radiale, donnent des termes oscillatoires indéfinis, qui sont traités à l'étape suivante. * Formation de l'expression finale pour l'équation matricielle d'élément. On substitue dans l'équation (30), qui est l'équation matricielle d'élément, toutes les expressions ci-dessus, comprenant à la fois les intégrales bien définies et les termes oscillatoires indéfinis. L'expres- sion entre parenthèses pour les trois matrices devient, [K] + io[C] o2[MW] = [K] - 2[] (32) + termes oscillatoires indéfinis et dans cette expression [K'] et M] J comprennent toutes les intégrales bien définies. Les termes oscillatoires indéfinis s'annulent tous, en laissant pour la forme finale de l'équation matricielle d'élément, [[K] P2[M]] (P} = (D} (33) avec les expressions suivantes pour K'. et M.i (les D. sont égales à zéro dans virtuellement toutes les applica- tions pratiques): K7 = B [A(}) R() + AM2 R(2) -2 A?) R(3)] = pr2 [A(,) R() - ú2 A(4) RJ (34) L vv 4' v1 p.. (4) avec les notations suivantes B est le module de compres- sibilité, ? est la densité, r1 est le rayon aplati du sphéroide d'élément infini et i1 (= f/r1) est l'excentri- cité du sphéroide d'élément infini. Les intégrales angulaires, ( i = 1,...,4, v',v = 1,...,n (voir le Tableau I), sont: A(}) = j'a asinJdd Vv v ' - -sin +9Jd d-q J> [a, ao ap a Jnr aVa ay av',Jav' 1 A, O0 IV sin O + Jddijn A vQ4) a sin Oaasin3 JOJdd i5 ao, Oq aóa wa, gsin3OJddl (35) avec = _ --|t dd pour des quadrilatères faO1fOddrI pour des triangles (36) Les quatre intégrales de l'équation (35) peuvent toutes être intégrées de façon numérique en utilisant des règles de Gauss classiques, du fait que les fonctions à intégrer n'ont pas de points singuliers et sont bornées(y compris le terme en 1/sinG, du fait que J-G9 lorsque 9 - 0). On calcule le jacobien, J, à partir de la projection de coordonnées de l'équation (7), ae au an1 aru (37) et on évalue les dérivées de ya, de la manière classique pour les éléments finis, a-DYV-a [] av; (38) L'angle G, pour les fonctions sin G et cos G, est également calculé à partir de l'équation (7). Les intégrales radiales R(',i=1.....4, ',=.....,m (voir le Tableau I), sont: R) = L^. (2Gm - 2 2ei2 + fi b2+2/ 3=2 2m RM3 La p2 - c:J+2 + R() 3=2 R(3) Rt'p = Lgt 12 cJ+2 B=2 R' = LW' ±F2H c2 =4 L1.4 r22- +C2e2mC+2JP(39) 2 L=1 avec L, = (1/k)ek(+" (40) = k-rl (41) = r e2 >2 1 (kr) -i2 o(x+)[t+ - Oe (x+ )-+E1 (cos 2x - i sin2x)dx (42) et ce sont des transformations de Fourier en sinus et en cosinus que l'on peut évaluer par des algorithmes classi- ques disponibles dans de nombreux progiciels mathémati- ques, _i2kr f J5 = kdr 2 1 1 - (kr) (43) r (43) e= f2 |cos 2x - sin 2x o 4(x+)2 _ Eú 2 (x+) - G = f - f-1 eJi2 ' kdr (44) ei2 |[(x+q)2 e2 2 - = e2 [ (x+)2 -1E (cos 2x - i si 2x) dx O+ i -1 H = |U e-i2k kdr (4 5 = e 2(o -) t x 2 2 [cos 2x -i sin 2x dxo O Elles peuvent toutes êtres évaluées avec le logiciel utilisé pour I Les relations entre ces intégrales sont: I3 = J3 -(kf)2 J+2 > - 1 (46) et H = G + (kf)2 J2 (47) De plus, aa = -ihga - (cc - 1)h_1 = 1,2....,m = 1,2,..., m+1 (48) O eth0 = hl.h+l = 0. = l =? i (49) et, d'après l'équation (48), aua = O pour Am + 1. C = ( hh, _ 7=1(0 et les h sont déterminés par l'équation (20). Les produits tensoriels Av'v RL'g dans l'équation (34) sont des matrices de N x N que l'on construit en multipliant chaque terme dans la matriceCR] de m x m, par la matrice [A], de n x n. entière, comme indiqué dans l'équation (51). Au-dessus et à gauche de la matrice se trouvent des tableaux, reproduits à partir du Tableau II, montrant la relation des index P et V avec j, et des index u' et la' avec i: 1...... t I 2À V 12... 1 2... I 2 V' i I' v' I I 2 2 R,,[A] R,2[A] R,M [A] 2 2 R:2[A] R,,2[A] ....DTD: ..................DTD: 1l5. v 'u^ 2 symétrique. R. [AI (51) L'équation (33), qui est l'équation matricielle d'élément, et les équations (34)-(51) sur lesquelles elle est basée, sont les équations mentionnées dans la case 20 de la figure lA. Ces équations sont suffisantes pour mettre en oeuvre l'invention. Elles peuvent être codées en n'importe quel langage de programmation scientifique; le FORTRAN est actuellement le choix le plus courant. On trouvera ci-après certaines observations pratiques concernant la résolution des équations ci-dessus par logiciel. Les intégrales radiales dans l'équation (39) sont identiques pour tous les éléments infinis dans une maille (du fait qu'elles sont indépendantes de variables angulaires et sont orientées le long de chemins radiaux aplatis identiques, c'est-à-dire des hyperboles co-focales partant du même sphéroïde aplati), ce qui fait qu'il suffit de les évaluer une seule fois pour un problème donné; leut coût de calcul est totalement négligeable. Par conséquent, l'intégration numérique qui est exigée pour générer CK'e] et LM-] pour chaque élément infini fait intervenir seulement l'évaluation des intégrales angulai- res 2-D, ce qui rend la génération de ces éléments 3-D aussi économique que celle d'éléments 2-D. De plus, du fait que la dépendance de l'élément vis-à-vis de la fréquence est contenue seulement dans les intégrales radiales, la génération d'éléments au cours d'un balayage de fréquence est pratiquement gratuite après la première fréquence. Il faut également noter que pour mettre en oeuvre la méthode des éléments infinis, il est nécessaire de générer unemaille 3-D dans le fluide entourant la structure. Ceci est un travail relativement simple pour des structures ayant des surfaces extérieures convexes. Deux ou trois ordres seulement seront exigés de façon caractéristique pour générer la maille de fluide entière pour une structure ayant une telle surface. L'utilisateur projettera simplement vers l'extérieur la maille 2-D qui se trouve sur la surface de la structure qui est en contact avec le fluide. Pour des configurations géométriques structurales plus générales, on peut générer des mailles 3-D par une grande variété de programmes de génération de maille qui sont commercialisés. L'étape de subdivision peut en outre comporter la construction d'une représentation géométrique du corps et la subdivision de la représentation du corps en éléments finis. En outre, le sphéroïde de limitation peut être un sphéroïde aplati minimal qui entoure la limite intérieure du fluide. Suivant un mode de réalisation préfère, l'étape de calcul de coefficients de matrice d'élément comprend l'application de l'équation de Helmboltz à une approximation d'un développement multipolaire d'une variable de champ acous- tique, cette approximation ayant la propriété de converger vers la solution exacte de l'équation au fur et à mesure que la maille est affinée. TABLEAU II Configuration de numérotation de noeuds liant l'index nodal (degré de liberté) j à l'index radial A et à l'index angulaire 0, N = m x n N=mxn j.1V J __ sur une surface sphéroidale 2 1 2 à r1 (base de l'élément infini) n i n n+l 2 1 sur une surface sphéroidale à r2 n+2 2 2 2n 2 n (m-!)n+ 1 m 1 (m-1)n+i m i sur une surface sphéroîdale à r (m-1)n+2 m 2 m N m n AN N EX E 1. Expressions formelles On calcule les matrices d'élément en partant d'un élément de taille finie, c'est-à-dire avec la face extérieure sur un sphéroïde de rayon aplati r, et en prenant ensuite la limite lorsque f o. Ainsi, en appliquant la méthode du résidu pondéré de Galerkin à l'équation (8) sur un seul élément, on obtient: Iim | | f (BV2p+co2pp)idV=O, i = 1,2,.....,N, (52) en utilisantk2= /BLa première intégrale est convertie en une intégrale de surface et en une autre intégrale de volume, en utilisant l'identité ViV2p='(ViVp)- Vyi'Vp et le théorème de la divergence. La substitution de l'équa- tion (11) dans les intégrales de volume donne l'équation matricielle d'élément suivante: ([K] - 2[M])(P} = (F{, (53) dans laquelle la matrice de raideur LK], la matrice de masse [M] et le vecteur de gradient de pression {F} sont respectivement, Kij = lim | f f BVViVVjadV, v(" M1i; = im f f f PXIjdV Fi = lim i Byi nap dS. (54) L'intégrale de surface-pour F est étendue sur la totalité de la limite de l'élément, S(e) Elle se décompose en deux intégrales: l'une sur la face extérieu- re S(e) et l'autre sur les faces restantes S(e) - s(e) Fi= lim | | Bv dS + lim Bvip dS. ion sa an; se -s" an (55) On considère la première intégrale dans l'équa- tion (55). Lorsque -a, les coordonnées sphéroidales aplaties tendent vers des coordonnées sphériques, ce qui fait que 3p/ Dn -- p/ Dr et dS -> r sin OdGdO, en désignant par r, 9, X des coordonnées sphériques. Pour évaluer Dp/ Dr lorsque f ->oo, on substitue les équations (11), (12) et (16) dans l'équation (10), ce qui donne une forme plus forte de la condition de Sommerfeld. aaP + ikp = o 1 r -, 0 (56) ar tr2J'(6 dans laquelle le O majuscule, qui se lit "grand O", signifie "d'ordre au moins égal à". On substitue %p/ Or dans l'équation (56) à)p/ cn dans la première intégrale. Le terme O(1/r) n'apporteaucune contribution à l'intégrale, du fait que Yi est O(1/r) et dS xr lorsque -Ao. Dans le terme ikp restant, on substitue l'équa- tion (11) pour le terme p. Dans la seconde intégrale dans l'équation (55), on substitue le reste de la quantité de mouvement linéaire à a p/ a n, c'est-à-dire Zp/ an = 2 [un, en désignant par u l'amplitude de la composante normale du déplacement n de la particule. En effectuant ces substitution, l'équation (55) devient: {F} = - ia)[C] (P + {D} (57) avec Cij = pclim f| yivjdS -- ± St,) Di = (Cpc)2 lim |f viudS (58) w X(4, (5) La matrice [C] est la matrice d'amortissement de rayonnement, représentant la perte de rayonnement jusqu'à l'infini. Le vecteur {D permet de spécifier un sur les faces latérales ou inférieure de l'élément. Ces faces sont soit sur lalimite du domaine, soit à l'intérieur du domaine. * Face sur la limite du domaine: une valeur différente de zéro correspond à une limite mobile avec une ampli- tude de déplacement (ou vitesse) normale connue. Le long des faces latérales, un terme u différent de zéro doit n décroître suffisamment rapidement suivant la direction radiale pour que l'intégrale dans l'équation (58) soit finie. Une valeur de zéro correspond à une limite rigide, c'est-à-dire qu'un mouvement est autorisé paral- lèlement à la limite, mais non perpendiculairement à celle-ci. Le cas le plus courant d'une limite infinie rigide est probablement celui dans lequel il y a un plan de symétrie à travers la structure et le fluide, permet- tant à l'analyste de modéliser seulement la structure et le fluide d'un côté du plan, c'est-à-dire que l'on a un problème de demi-espace infini. Un plan de symétrie apparaît lorsque toutes les propriétés passives du système (configuration géométrique, contraintes, propriétés des matériaux et, éventuellement, précon- traintes hydrostatiques) sont symétriques par rapport à un plan. Si l'excitation est également symétrique, ou si on considère seulement la composante symétrique d'une excitation asymétrique, alors la réponse du système doit posséder la même symétrie. Par conséquent, des quantités vectorielles peuvent être symétriques par rapport au plan. Dans le cas du déplacement de particules, ceci signifie que la composante normale doit s'annuler, afin de préserver la continuité de la matière. * Face à l'intérieur du domaine: Toute face intérieure est commune à deux éléments acoustiques infinis ou de taille finie adjacents. L'assemblage des vecteurs [D5 provenant de deux éléments adjacents donne, sur l'inter- face, une somme de deux intégrales de surface avec un signe opposé, du fait que les normales dirigées vers l'extérieur sont dans des directions opposées. La paire d'intégrales représente donc la différence de u de part n et d'autre de l'interface. Une différence non nulle correspond à une discontinuité sous la forme d'un saut dans le déplacement normal, c'est-à-dire une séparation du fluide. Une différence nulle correspond à la préser- vation de la continuité de la matière, ce qui est la condition habituelle dans des applications pratiques. En substituant l'équation (57) dans l'équation (53), on obtient pour l'équation matricielle d'élément: ([K] + iW[C] - 02 [M])(P) = (D) (59) dans laquelle, pour résumer, les matrices de raideur, de masse et d'amortissement de rayonnement sont respective- ment, Kij = lim S BVvi ÀV jdV r--±V( Mij = lirm f JJJ pivjdV V(,) Ci; = pclim f fi4jdS (60) Dans ce qui suit, on tranforme les intégrales dans l'équation (60) en expressions qui peuvent être évaluées de façon numérique. Le vecteur D} est égal à zéro dans presque toutes les applications pratiques. 2. Transformation d'intégrales en coordonnées sphéroîdales aplaties Les éléments différentiels de volume et de surface sont dV = J vdrd9do et dS = JsdGdo, en désignant par JV le jacobien de volume et par JS le jacobien de surface (sur des sphéroïdes). Du fait que les faces de l'élément infini se conforment à des sphéroïdes (surfaces G, 0) et à des hyperboles (lignes r), les limites d'inté- gration pour les intégrales de volume peuvent être séparées en limites angulaires et en limites radiales. L'équation (60) devient donc: Kij = lim f f f BVyi VjJvdrdedo r-*O) r1 Mij = lim f pfifjJdrdedO Rico l'élément, c'est-à-dire toute surface sphéroidale co- focale intérieure à l'élément et limitée par les faces latérales. Il n'y a qu'une seule section sphéroîdale dans l'espace r, 9, 0, du fait que les coordonnées 9, 0 de la limite de la section transversale sont indépendantes de r. Par consé- -(e) éuvuà (e) quent, S(e) équivaut à (e) dans l'intégration de surface pour Cij. La transformation de JV' JS et V en coordonnées sphéroidales allongées est la suivante. * Jacobien de volume ad ad ad JV = a X = hra, 'hoaox hao = hheho (62) ar ae -( avec les notations suivantes: d = xi + yj + zk est le vecteur de position et ar, a9, a sont des vecteurs uni- taires, avec des dimensions de longueur, respectivement suivant les directions r, G, 0. Pour des coordonnées orthogo- nales, on a ar a x a0 = 1. Les facteurs d'échelle hr, hG et h sont donnés par: ad ax ay az h =- =- /-a + + u= r,O, ( 0-u X au au Dau ' (63) L'utilisation de l'équation (5) donne: hr / r _ f2 sin20 r2 - f2 hE = fr 2 f2 sin20 ho = r sin 0 (64) et l'équation (64) devient r2 _ f2 sin2 O -v r sin (65) /r2 _f2 * Jacobien de surface J est nécessaire dans Cii seulement sur la S 13 surface extérieure, c'est-à-dire sur le sphéroïde à r = Js(r) = ad ad haa x ha0a = he(P)ho( ) I = r2 _f2 sin20 sinO (66) et il est nécessaire seulement à-la limite lorsque - (lorsque f reste constant). Par conséquent, lim Js = 2 sinO (67) i --.} ce qui est identique au jacobien de surface en coordonnées sphériques, du fait que les coordonnées sphéroîdales aplaties tendent vers les coordonnées sphériques lorsque ? --. * Gradient h, ar + ao ho a + ao ho (68) Le gradient est nécessaire seulement dans K., o il i] apparaît sous la forme Vli' VyjJv. En utilisant les équa- tions (64), (65) et (68), on obtient: Vyi ' VIijv = r[r2 - f2 r sinOasinav ar or r r 2 a_ f2 ae + r2 _ f2sin2O a1i aV1j r r2 _ f2 sinOa aO (69) On peut maintenant substituer les équations (69), (65) et (67) respectivement dans Ki.., Mij.. et Cij, dans l'équation (61). De plus, on substitue l'équation (12) à Yi et '4fj, en introduisant le symbole prime pour distinguer les index V et p qui sont associés à l'index i: 7i(,11,r) = va,(,r)Tt,(r) ici(,rj,r) = (,r)7t(r) (70) En effectuant ces substitutions, et en notant que les fonctions de forme angulaire sont implicitement des fonctions de G, 0 (à cause de la projection de coor- données dans l'équation (7)), toutes les intégrales 3-D se (i) i séparent en produits d'intégrales "angulaires" 2-D, A(,i), et d'intégrales "radiales" 1-D, R(,i. L'intégrale de PP surface 2-D pour Cij se sépare d'une manière similaire. Les expressions résultantes sont les suivantes: K= Br[A)O(1) + - _ 2 AV2 R Mij = BpAA2 p('(4) _ ú2 A() R(3J Mi Pr[A v(v)Rg'g + vv,.. vvVg V'V r,. Il. 71 Cij =-pcA() lim ur r(7 avec A(,A[ = ff, vva 'sinedOdo (i(') =|| av av sin a a8/ i aoa, =: a( av in sn sO + i d a av A v3) = y' be ayv sin Odedo V' A v?,) = || rvsin3'OdOd (72) ! I, im r r2 _f2 d' dwg --dr -m= r2 dr dr rr r R2) = lim óf'xdr - -oo -'4 r2 -f2 R(4) = lim r r r2 -f r (73dr r. /2 _f2 =i.r2 r 1 - (73) et - = excentricité du sphéroïde d'élément infini 1r (74) Les symboles R(1), et R(4) sont surmontés par une barre du fait qu'un terme sera ultérieurement éliminé de chaque intégrale (voir les équations (87) et (101)), et que les termes restants seront désignés par Rlu) et R() afin que les quatre intégrales radiales aient toutes des symboles uniformes (sans barres supérieures) dans les expressions finales (équation (112)). Les produits tensoriels A., R, sont des matrices N x N qui sont construites en multipliant chaque terme dans la matrice [RJ m x m, par la matrice [A], n x n, entière, comme représenté dans l'équation (75). Au-dessus et à gauche de la matrice se trouvent des tableaux, reproduits à partir du Tableau II, montrant la relation des index p et V avec j. et des index j' et V' avec i. j 1 2. ..... v A| I 2.. V 1 2 n 1 2.n 1 2 n i ' v' 2 2 i R11[A] R12[A] RIm[A] l........................... ............................................ 2 2......... R22[A]R2m, [A] symétrique Rmm, [A] (75) * Intégrales angulaires Du fait de la projection de coordonnées de l'équation (7), on a dGdo = Jd d 1, avec pour le jacobien: ae au J=a a3 ae a3 auxq ar(76) Par conséquent, l'équation (-72) devient, ail A(1 fi |V':v sinO Jd{ dl A. V- sin o + au sO 04a) 30 3ó sin O A( f si= nm nJdn dri A^) = S| Va VIsin3'Jdé du..DTD: "(77) avec o(') d {''' dé dT pour des quadrilatères =f) f1l.o!.. dér pour des triangles (78) Les quatre intégrales peuvent toutes être intégrées de façon numérique en utilisant des règles de Gauss, du fait que les fonctions à intégrer n'ont pas de points singuliers et sont bornées (y compris le terme 1/sinG, du fait que J oc O lorsque G -- 0). On évalue les dérivées de 8(a) de la manière classique, ao = ac J l@ J (79) En fait, la totalité de l'évaluation des intégrales angu- laires suit les procédures classiques d'analyse par la méthode des éléments finis. 4. Intégrales radiales -(1) * " D'après l'équation (16), d _ -i-k(r-r) m+ aa = ke- ik(r-rJ Y. ga(80) dr cc.=1(kr)8 avec aga = -ih-a - (l-)ha. = 1,2...m = 1,2.m +1 (81) et h.o = hg,m+ l = O. di4tf dr k2eL k(r, + r) [2m+2 b5 ei2k dr dr 5 =2 (kr) (82) avec *-1= (83) -Y= (83) et d'après l'équation (81), on a a = O pour ">m + 1. En substituant l'équation (82) dans la première des équations (73), on obtient: 2m f - R(,) = La im 2, lf iraed + 2m I e-i2krk dr (kr) e (84) avec Llk (85) Un facteur k est introduit dans les intégrales (et est annulé par le 1/k dans L 1) pour les rendre sans dimen- sion. L'intégrale b2 est séparée des autres du fait que la limite supérieure donne un terme indéfini (qui sera annulé ultérieurement par des termes similaires provenant de M. et C..ij). D'une manière formelle, on écrira: lim |, r -1e kdr = r - f1ei2"k d ire - i2; i -i2kr, + hm-e-e-- __-..2 2 (86) Les autres intégrales sont bien définies, ce qui fait que le processus de limitation exige seulement de changer la limite supérieure en la faisant passer de f à o. L'équa- tion (84) devient (1)-L b2iUm e + R() (87) avec R( L 2G e-i2 + b L j=l (88) (89] = krl, (89) et 1 G [ - i e i2kr k dr, (90) f (kr)g = (k' e-L2A2 21 (91) On peut évaluer les intégrales G et I1 en faisant le changement de variable: x = kr - (92) Les quations (90) et (91) deviennent(92) Les équations (90) et (91) deviennent G=e'24 fe { (x+4)2- _ j (cos2x - i sin 2x)dx3 (9 3 - E2 e-i24 ( +42 _e2 o = e'I (x4+l (cos2x - i sin 2x)dx (94) qui sont des transformées de Fourier en sinus et en cosinus que l'on peut évaluer par des algorithmes standards disponibles dans de nombreux progiciels mathéma- tiques. * R()m D'après l'équation (16), It t = LW Y. ei2 V L'WL= Lej ( [) k (95) avec -! CB = z hw'y(96) Y= I et h = O pouro > m. En substituant l'équation (95) dans la seconde des équations (73), et en prenant l'infini pour la limite supérieure, du fait que toutes les intégrales sont bien définies, on obtient 2m R(2,)= Lg c (97) =2 273437g avec je e- i[e r kdr 3> 1 (Jc)r (98) Le terme 1 - (f/r) est supérieur à 0, du fait que r r1 > f. L'utilisation de l'équation (92) convertit J en transformées de Fourier en sinus et en cosinus, JF = e-i2 È cos2x- isin 2x d 0 (x+;)2 _ s2;2 (x+;)-l (99) qui peuvent être évaluées pa-r le progiciel mathématique qui est utilisé pour évaluer G et I A ci-dessus. * R(3)m En substituant l'équation (95) dans la troisième des équations (73) et en prenant l'infini pour la limite supérieure, du fait que toutes les intégrales sont bien définies, on obtient R(3) 2 R,?= Lp., E CPJ+2 (100) :=2 f 4* M M En substituant l'équation (95) dans la quatrième des équations (73), et en procédant de la même manière que pour R(1) (c'est-à-dire en séparant l'intégrale C2 des autres, du fait que la limite supérieure produit un terme oscillatoire indéfini), on obtient -(4) 1 C2i R = L 2 li2 m e2 + R(4)101) 2 (101) avec c2 i' 7;{2m -2 R:)L t H -C21 i2ç + E2z =L'j F2 H -+ c+ 2J1 49 2 j=1 (102) et i - (103) Hf 2 ei2r k dr(13 Le fait d'utiliser à nouveau l'équation (92) convertit H en transformées de Fourier en sinus et en cosinus, H=e- i2j [ - 1 (cos 2x-i sin 2x)dx o (x+)2 _ l2 (104) que l'on peut également évaluer par le progiciel mathéma- tique utilisé pour évaluer G, I et J, ci-dessus. 5. Fonctions radiales sur la face extérieure La matrice d'amortissement de rayonnement, Ci., dans l'équation (71) contient le terme lim g,' r En utilisant l'équation (95), on obtient, [n = n k 2, =2 C5 rs [ j 2r k- 5L=2 (e) 2 (105) Le fait de prendre-la limite lorsque f -y X élimine tous les termes co, à l'exception de c2, qui contient un terme oscillatoire indéfini: lim [1rv'J;2 = L"'>.T i. e- l. [lv,] î2r (02 lim)e-i2) - (106) 6. Disparition de termes oscillatoires indéfinis On peut maintenant substituer dans les matrices de raideur, de masse et d'amortissement, dans l'équation (71), les expressions pour toutes les intégrales radiales (équations (87), (88), (97), (100), (101) et (102)) et pour les fonctions radiales sur la face extérieure (équation 106). En additionnant ensemble ces trois matri- ces, comme dans l'équation matricielle d'élément, c'est-à- dire l'équation (59), on obtient [K] + i w[Ci - w [M] = [K] - u 2[M + termes oscillatoires indéfinis (107) Dans l'équation ci- dessus, [K' ] et [M] comprennent toutes les intégrales bien définies et sont résumées au paragraphe suivant. Les "termes oscillatoires indéfinis" comprennent tous les termes qui contiennent l'expression - i2kf lim e Ces termes apparaissent en trois occasions, correspondant aux équations (87), (101) et (106), qui contribuent respectivement à [K], [C] et [M]. On a donc, termes oscillatoires indéfinis = BAVL 2 lirm ei2 + icopcA(l)L, li -i2ke BCpA(vl') -W 2Iim o2pr2A()L ci i -i2k^ - I V-lime A'v I- g' 22;--,e = (b2 + C2)B-2 A(\) L lim e-i2k; (b 2 +; v v Lg,g --*- (108) en utilisant c = ck et c = B/. D'après les équations (81) et (83), 2 = ap lal = - h,,lh1, et d'après l'équa- tion (96), c2 = hP, lh = - b2. Par conséquent, termes oscillatoires indéfinis = 0 (109) Sous la forme de phrases, ceci s'exprime de la manière suivante: les matrices de raideur et de masse consistent en deux types de termes. Ceux qui sont indépen- dants de -la position de la face extérieure (les intégrales "bien définies") et ceux qui dépendent effectivement de la position et qui oscillent donc lorsque la face s'éloigne jusqu'à l'infini. La matrice d'amortissement, qui repré- sente l'application de la condition de rayonnement de Sommerfeld à la face extérieure, est complètement oscil- latoire. L'équation (109) indique que la condition de rayonnement annule exactement les termes oscillatoires dans les matrices de raideur et de masse. 7. Expressions finales Il résulte des équations (107) et (109) que la forme finale de l'équation matricielle d'élément pour l'élément acoustique infini est [K]- W2[Mjl[P = {D) (110) avec K" = B A(+) R1) + A(2) R) - ú2A(3) R(3) M-j = pr? [A$l)R:) - 2(4, (111) M [ =( pr [A) R( A, o2) R(,] Le vecteur {D} est donné par l'équation (58); il est égal à zéro dans presque toutes les applications pratiques. Les intégrales angulaires A(i) i = 1, 2, 3, 4, sont données VIV par l'équation (77). Les intégrales radiales sont R1) 2 La A [L2G- Lb2&L2 m +m b+. W'4 g' 2 2-b? + b+2IB R = L" 2m =2 2m R(3) -L2 Rd = L4,g4 Y C5J++2 P=2 R = Lpt 4 C2H - 'c2e-i2 m =1 (112) et L, G, I J, b, c3, et sont donnés respectivement par les équations (85), (90), (91), (98), (103), (83), (96) et (89). Les équations (83) et (96) nécessitent à leur tour les h, qui sont donnés par l'équation (20). Les intégrales radiales dans l'équation (112) sont identiques pour chaque élément infini dans une maille (du fait qu'elles sont indépendantes de variables angu- laires et sont orientées le long de chemins radiaux allongés identiques, c'est-à-dire des hyperboles qui partent du même sphéroïde), ce qui fait qu'il suffit de les évaluer une seule fois pour un problème donné; leur coût de calcul est totalement négligeable. Par conséquent, l'intégration numérique qui est exigée pour générer [K] et [M'"] pour chaque élément infini nécessite seulement l'évaluation d'intégrales angulaires 2-D, ce qui rend la génération de ces éléments 3-D aussi économique que celles d'éléments 2-D. De plus, du fait que la dépendance de l'élément vis-à-vis de la fréquence est contenue seulement 273437g dans les intégrales radiales, la génération d'éléments pendant un balayage de fréquence est pratiquement gratuite après la première fréquence. REVENDICATIONS 1. Procédé pour faire fonctionner un ordinateur, ayant au moins une mémoire numérique et au moins un élément de traitement de données, pour simuler le compor- tement acoustique d'un corps entouré par un milieu fluide, ce corps ayant une surface extérieure et étant soumis à des conditions d'excitation données, caractérisé en ce qu'il comprend les étapes suivantes: (a) on subdivise au moins le milieu fluide en une configuration d'éléments, cette configuration étant appelée une maille, et on enregistre cette maille dans la mémoire; (b) pour chaque élément de la maille, on calcule un ensemble de coeffi- cients de matrice d'élément, et on enregistre ces coeffi- cients dans la mémoire; (c) on assemble tous les coeffi- cients de matrice d'élément en une matrice de système, et on enregistre la matrice de système dans la mémoire; (d) dans l'élément de traitement de données, on résout l'équa- tion matricielle de système, pour créer ainsi une descrip- tion des valeurs prises par une variable de champ acous- tique dans toute l'étendue de la maille; et (e) on enre- gistre la description dans un dispositif d'enregistrement de données, dans lequel l'étape de subdivision comprend les étapes suivantes: (f) on construit une limite intérieure du milieu fluide, de façon que cette limite intérieure coïncide avec la surface extérieure du corps; (g) on construit un sphéroïde aplati, que l'on appelle le sphéroïde de limitation, qui enferme la limite intérieure du fluide; et (h) on remplit un espace entou- rant le sphéroïde de limitation avec des éléments, que l'on appelle des éléments infinis, avec les conditions suivantes: (i) chaque élément infini est limité par une base, au moins trois faces latérales et une face extérieu- re; (ii) chaque base respective s'étend sur le sphéroïde de limitation; (iii) chaque face extérieure respective appartient à une surface sphéroldale aplatie qui est co- 273437g focale au sphéroïde de limitation; (iv) chaque face latérale respective est un lieu d'hyperboles co- focales au sphéroïde de limitation; et (v) la face extérieure s'éloigne jusqu'à un rayon aplati infini. 2. Procédé selon la revendication 1, caractérisé en ce que l'étape de subdivision comprend en outre la construction d'au moins une couche intermédiaire d'élé- ments finis entre la limite intérieure du fluide et le sphéroïde de limitation. 3. Procédé selon la revendication 1, caractérisé en ce que le sphéroïde de limitation coïncide avec la limite intérieure du fluide. 4. Procédé selon la revendication 1, caractérisé en ce que l'étape de subdivision comprend en outre: (a) la construction d'une représentation géométrique du corps; et (b) la subdivision de la représentation du corps en éléments finis. 5. Procédé selon la revendication 1, caractérisé en ce que le sphéroïde de limitation est un sphéroïde aplati minimal entourant la limite intérieure du fluide. 6. Procédé selon la revendication 1, caractérisé en ce que l'étape de calcul de coefficients de matrice d'élément est accomplie de manière à satisfaire la condi- tion de rayonnement de Sommerfeld. 7. Procédé selon la revendication 1, caractérisé en ce que l'étape de calcul de coefficients de matrice d'élément comprend l'application de l'équation de Helmholtz à une approximation d'un développement multi- polaire d'une variable de champ acoustique, cette appro- ximation ayant la propriété de converger vers la solution exacte de l'équation au fur et à mesure que la maille est affinée. 8. Machine pour simuler le comportement acous- tique d'un corps entouré par un milieufluide, le corps ayant une surface extérieure, et étant soumis à des conditions d'excitation données. caractérisé en ce au'elle ccrmnrend: (a) des moyens pour subdiviser au moins le milieu fluide en une configuration d'éléments, cette configuration étant appelée une maille; (b) un élément de mémoire numérique pour enregistrer la maille; (c) des moyens de traitement numérique pour calculer un ensemble de coefficients de matrice d'élément pour chaque élément de la maille; (d) un élément de mémoire numérique pour enregistrer les coeffi- cients de matrice d'élément, assemblés à partir de tous les éléments, sous la forme d'une matrice de système; (e) des moyens de traitement numérique pour résoudre la matrice de système, pour créer ainsi une description des valeurs prises par une variable de champ acoustique sur toute l'étendue de la maille; et (f) des moyens pour enregistrer la description résultante de la variable de champ acoustique, dans laauelle les movens de subdivision comprennent: (g) des moyens pour construire une limite intériere du milieu fluide, de façon que cette limite intérieure coïncide avec la surface extérieure du corps; (h) des moyens pour construire un sphéroïde aplati, que l'on appelle le sphéroïde de limitation, qui entoure la limite intérieure du fluide; et (i) des moyens pour remplir un espace entourant le sphéroïde de limitation avec des éléments, que l'on appelle des éléments infinis, avec les conditions suivantes: (i) chaque élément infini est limité par une base, au moins trois faces latérales et une face extérieure; (ii) chaque base respective s'étend sur le sphéroïde de limitation; (iii) chaque face exté- rieure respective appartient à une surface sphéroidale aplatie co- focale au sphéroïde de limitation; (iv) chaque face latérale respective est un lieu d'hyperboles co-focales au sphéroïde de limitation; et (v) la face extérieure s'éloigne jusqu'à un rayon aplati infini. 9. Machine selon la revendication 8, caracté- risée en ce que les moyens de subdivision comprennent en 273437g outre des moyens pour construire au moins une couche intermédiaire d'éléments finis, entre la limite intérieure du fluide et le sphéroïde de limitation. 10. Machine selon la revendication 8, caracté- risée en ce que le sphéroïde de limitation coïncide avec la limite intérieure du fluide. 11. Machine selon la revendication 8, caracté- risée en ce que les moyens de subdivision comprennent en outre: (a) des moyens pour construire une représentation géométrique du corps; et (b) des moyens pour subdiviser la représentation du corps en éléments finis. 12. Machine selon la revendication 8, caracté- risée en ce que le sphéroïde de limitation est un sphéroïde aplati minimal entourant la limite intérieure du fluide. 13. Machine selon la revendication 8, caracté- risée en ce que les moyens destinés à calculer des coeffi- cients de matrice d'élément sont soumis à des contraintes de façon à satisfaire la condition de rayonnement de Sommerfeld. 14. Machine selon la revendication 8, caracté- risée en ce que les moyens destinés à calculer des coeffi- cients de matrice d'élément comprennent des moyens destinés à appliquer l'équation de Helmholtz à une appro- ximation d'un développement multipolaire d'une variable de champ acoustique, et en ce que cette approximation a la propriété de converger vers la solution exacte de l'équa- tion au fur et à mesure que la maille est affinée.