Modélisation du dégivrage aéronautique
Modélisation du dégivrage aéronautique
Présentée par
Je voudrais commencer par mes trois encadrants de thèse, Edith, Nicolas et Fred,
merci pour les années passées à échanger sur ces recherches et pour l’opportunité et la
confiance que vous m’avez donnée en m’acceptant dans votre équipe. Je garderais une
grande leçon et le souvenir des moments positifs et agréables du projet. Merci à vous
d’avoir partagé vos connaissances avec moi, pour avoir fait preuve d’engagement et
patiente jusqu’à la dernière phase de la thèse.
Je tiens à remercier M. Olivier Chadebec pour votre soutien et conseil lors des CSI
de fin de ma thèse, avec ceux de Yves Lembeye, vos remarques et propositions si
pertinentes ont été un élément clé pour l’aboutissement de mes travaux de rédaction.
À tous les membres de l’équipe EP, enseignants, doctorants, stagiaires, post doc
des différentes générations que j’ai pu croiser : merci pour ce partage ! j’ai toujours passé
des super moments dans la bonne ambiance qui règne chez nous, là où ça sent le
composant cramé oui.
Merci à l’ensemble du laboratoire G2Elab, à l’équipe informatique et technique, mes
copines du service RH et du service financier (petite pensée pour Florence, Nadine, Cathy
et pour toi ma Lili <3), merci pour votre écoute et votre travail dévoues.
Un souvenir spécial pour les grandes amitiés que j’ai pu cultiver au long de ces
années au labo : Manel, Geneviève, Johan, Fabien, Gatien, Guillaume, Ekaitz, Antoine,
Selle, Stéphane, Jeremy, Victor, Yafang, Hakim, Max, Clem, Guillerme, Antho, Bob,
Vincent, Vihn, mon cousin Farshid, Yunxia, Nils, Sandrine, Olivier, Thibaut, Long, Becka,
Badr ma chérie, Andrés, Clara, Olga, Zaki, Baoling….. ufff vous êtes nombreux, j’ai de la
chance ! vous êtes le plus grand trésor que je garderais de cette expérience de 5 ans.
Mais c’est avant tout grâce à Dieux/l’Univers qui m’a concédé la santé, la force, la
musique, la paix et l’occasion de vous rencontrer tous, que je peux maintenant dire avec
joie et fierté que cet objectif est accompli. Une nouvelle étape se profile dans le
paysage…merci pour le souffle de vie.
Table des matières
IV.3.e Analyse des pertes par courants de Foucault parasites au niveau de la surface
à dégivrer ................................................................................................................ 111
IV.5 Assemblage des modèles et approche statique avec boucle de calcul................ 114
IV.6 Conclusions sur la modélisation électrique .......................................................... 116
V.2.b Etude de sensibilité avec la définition des éléments plots ......................... 131
V.3.b Essai expérimental 2 : le déplacement est imposé par les forces générées au
niveau des lignes des conducteurs ......................................................................... 139
Tableau 1-1. Classement de taille de cabine d’un avion selon la FAA [9]. ........................ 14
Tableau 1-2. Comparatif des technologies existantes. ..................................................... 23
Tableau 1-3. Descriptif de la constitution du tapis par couches matérielles [33]. .............. 30
Tableau 1-4. Propriétés des matériaux de construction pour le démonstrateur [33]. ........ 32
Tableau 1-5. Cahier de charges. ...................................................................................... 37
Tableau 2-1. Valeurs d’entrée pour les paramètres géométriques et physiques des
premières simulations pour la force.................................................................................. 51
Tableau 2-2. Valeurs d’entrée pour les paramètres géométriques et physiques du 1er
scenario. .......................................................................................................................... 58
Tableau 2-3. Valeurs d’entrée pour les paramètres géométriques et physiques du 2ème
scenario. .......................................................................................................................... 59
Tableau 2-4. Valeurs d’entrée pour les paramètres géométriques et physiques du 3ème
scenario. .......................................................................................................................... 61
Tableau 3-1. Traduction des conditions aux limites au niveau des grandeurs mécaniques.
........................................................................................................................................ 75
Tableau 3-2. Calcul comparatif de différents ∙ pour un nombre d’éléments n. ........... 79
Tableau 3-3. Description des trois cas comparatifs pour l’étude du feuilletage des
conducteurs. .................................................................................................................... 80
Tableau 3-4. Données géométriques pour l’analyse de l’impact du feuilletage des
conducteurs. .................................................................................................................... 81
Tableau 3-5. Cas étudiés pour la décomposition des conducteurs en multiples lamelles. 82
Tableau 3-6. Conditions de continuité pour une arête = . ......................................... 90
Tableau 3-7. Conditions aux limites au niveau des arêtes = et = . ..................... 91
Tableau 4-1. Bilan des énergies dans le système. ......................................................... 109
Tableau 4-2. Pertes totales dans la plaque dues aux courants induits liées à l’impulsion de
courant........................................................................................................................... 112
Tableau 4-3. Energie et pertes au niveau de la surface à dégivrer. ................................ 113
Tableau 4-4. Résume sur l’analyse énergétique de l’ensemble du système................... 113
Tableau 5-1. Comparaison entre la configuration initiale et la structure avec 6 conducteurs.
...................................................................................................................................... 124
Tableau 5-2. Comparaison entre le cas initial et la structure avec conducteurs de largeur
double. ........................................................................................................................... 126
Tableau 5-3. Comparaison entre le cas initial et la structure avec moins de plots. ......... 127
Tableau 5-4. Comparaison entre le cas initial et la structure avec un tapis moins large. 128
Tableau 5-5. Comparaison entre le cas initial et la structure avec des conducteurs plus
longs. ............................................................................................................................. 129
Tableau 5-6. Récapitulatif des résultats obtenus avec le modèle statique. ................... 1310
Notations
Modèle électromagnétique
𝐵 Induction magnétique [T]
𝐼 Courant électrique [A]
𝐹 Force de Laplace [N]
Modèle mécanique
𝛺 Domaine d’un milieu continu
𝑆 Déformation
𝜎 Contrainte
𝑢, 𝑤 Déplacement [m]
𝐶 Tenseur d’élasticité
𝜆∗ , 𝜇∗ Coefficients de Lamé
𝑅 Réaction
𝑁 Effort normal
𝑉 Effort tranchant
𝑀 Moment de torsion ou de flexion [°/m]
𝜌 Masse volumique [kg/m]
𝑝⃗ Force appliqué par unité de longueur [N/m]
Γ⃗ Couple par unité de longueur [N·m]
𝐼𝑧 Moment quadratique [m4]
𝜃 Angle de rotation [°]
𝐺𝑥𝑦 Coefficient de cisaillement
𝑊 Déformation d’une plaque [m]
Modèle électrique
𝑊, 𝐸 Énergie [J]
𝑉 Potentiel électrique [V]
𝐼 Courant électrique [A]
𝑃 Puissance [W]
𝐶 Capacité [F]
𝑅 Résistance [Ohm]
𝐿 Inductance [H]
𝑍 Impedance [Ω]
𝑇 Période [s-1]
Paramètres physiques
𝜇0 Constante de perméabilité magnétique dans le vide 4πx10-7 (H/m)
𝐸 Module d’Young
𝜐 Coefficient de Poisson
Introduction :
Enjeux de l’avion plus électrique
Introduction générale : Enjeux de l’avion plus électrique
La conception d’un avion plus électrique entraine l’élimination d’une partie importante des
tuyaux liés aux systèmes mécaniques, ce qui conduit à une réduction de poids et volume
et ainsi que de la consommation énergétique [2]. En contrepartie, remplacer les dispositifs
classiques par des dispositifs électriques implique la gestion des problèmes liés à la CEM,
ainsi que la recherche de solutions de stockage d’énergie et le développement de
convertisseurs beaucoup plus performants, nécessaires pour répondre aux nouveaux
besoins énergétiques.
3
Introduction générale : Enjeux de l’avion plus électrique
pratiquement éliminés, diminution du taux des victimes par passager par km d’un facteur
10 [3].
Projet CORAC-GENOME
La France, historiquement acteur clé de la recherche aéronautique, a lancé en 2013 le
projet Gestion Optimisée de l’Energie «GENOME». GENOME est un projet national de
grande envergure développé au sein du Conseil pour la Recherche Aéronautique Civile
«CORAC», organisme qui réunit les principaux acteurs de l’industrie aéronautique
française [4].
4
Introduction générale : Enjeux de l’avion plus électrique
a) b)
Figure ii - Givre blanc a) et verglas b) [6].
Pour se prémunir de ces effets, il existe des dispositifs de dégivrage au sol et des
dispositifs de protection qui opèrent pendant le vol. Les premiers consistent généralement
en un bain de substances chimiques qui altèrent le point de solidification de l’eau,
notamment le propylène glycol, avec une action qui dure encore quelques minutes après
le décollage, le temps que les systèmes d’antigivrage et de dégivrage en vol soient actifs.
5
Introduction générale : Enjeux de l’avion plus électrique
- modéliser le système et comparer les résultats aux tests sur les prototypes,
- à partir du modèle validé, identifier les paramètres de conception,
- trouver des pistes pour optimiser le dispositif (explorer les degrés de libertés
géométriques, proposer des stratégies de commande pour la source
d’alimentation, adaptation du modèle pour une future optimisation).
- l’électromagnétisme,
- la mécanique,
- l’électronique de puissance.
Pendant les trois années de la thèse, le travail développé s’est principalement concentré
sur la compréhension approfondie du phénomène de déformation élastique par action des
efforts électromagnétiques afin d’aboutir à un modèle multiphysique adapté permettant
d’identifier des paramètres d’optimisation.
6
Introduction générale : Enjeux de l’avion plus électrique
complexité du problème, d’enrichir les modèles dans un second temps selon les
précisions des résultats obtenus et le développement de l’analyse.
Les phénomènes physiques ont donc été étudiés séparément, créant des modèles
indépendants. Une deuxième étape a porté sur l’assemblage de ces blocs comme
éléments constitutifs d’une chaine de calcul. L’assemblage, achevé dans un premier
temps en statique, a permis ainsi de valider la cohérence de leur couplage. L’aspect
dynamique ajouté dans un second temps, a mis en évidence la difficulté liée à cette
modélisation multiphysique très contraignante.
- modélisation électromagnétique,
- modélisation mécanique,
- modélisation électrique.
Structure du manuscrit
Le manuscrit de thèse se divise en cinq chapitres :
7
Introduction générale : Enjeux de l’avion plus électrique
8
Chapitre 1 :
Le dégivrage en aéronautique et
ses technologies
Chapitre I : Le dégivrage en aéronautique et ses technologies
I.1 Introduction
Dans la suite, nous allons présenter un résumé de ces solutions pour essayer de mieux
comprendre leur mode de fonctionnement. Les industriels sont conscients que la solution
optimale pourrait se trouver dans un système hybride combinant plusieurs d’entre elles. Le
choix de la solution sera ainsi un compromis entre les différents critères : technologiques,
disponibilité des matières premières, coût, impact écologique, modularité et maintenance,
sécurité, robustesse, entre autres.
Cette solution consiste à traiter les surfaces susceptibles d’être givrées avec des
substances qui empêchent la glace de s’y fixer, appliqués au chaud. Il s’agit souvent de
fluides à base de glycol visant à se mélanger avec les particules d’eau en contact avec la
surface et modifier leur point de fusion vers une température beaucoup plus basse [1].
Usuellement on utilise des mécanismes externes à l’avion pour appliquer ces substances à
pression sur l’ensemble des surfaces de l’avion, mais principalement au niveau des bords
d’attaque des ailes, la dérive et aux rebords des turbines comme illustré sur la Figure 1-1.
Ces fluides contiennent des glycols, de l’eau et des substances anticorrosion. Par
ailleurs, pour le cas des fluides antigivrants, ils peuvent être mélangés avec des polymères
12
Chapitre I : Le dégivrage en aéronautique et ses technologies
pour leur conférer une meilleure viscosité ; ce qui se traduit par une meilleure tenue sur les
surfaces malgré les changements de pression.
Des normes ISO 11078 :2007 et AMS1428G définissent les standards actuels de
sécurité, qualité et emploi pour ces produits [4] [5]. Ces solutions sont combinées aux autres
méthodes pour accroitre leur efficacité. Il n’en demeure pas moins important de prendre en
considération la pollution liée à la synthèse et l’utilisation de telles substances.
Figure 1-2. Boudins pneumatiques dégivrants en caoutchouc sur le bord d’attaque d’une aile [6].
En présence de conditions givrantes, les boudins sont gonflés pendant quelques instants
pour expulser la couche de glace indésirable à travers le choc mécanique généré au
moment de son expansion. La Figure 1-3 représente les deux positions possibles.
1En aéronautique, un fluide dit Non-newtonien c’est quand il est mélangé avec des polymères pour
altérer ça viscosité, et Newtonien quand il ne contient pas des polymères ajoutés.
13
Chapitre I : Le dégivrage en aéronautique et ses technologies
a) b)
Figure 1-3 Boudins pneumatiques en néoprène en position initiale-dégonflés a) et actionnés b) [7].
Cette technologie est une des premières solutions contre le givre imaginées et elle est
utilisée encore de nos jours [8]. Cependant, cette solution est plus adaptée aux avions
turbopropulsés de taille petite à moyenne à cause des contraintes volumiques imposés par
le système pneumatique. On peut consulter le Tableau 1-1 pour voir le classement par taille
de cabine selon l’Administration Fédérale d’Aviation américaine.
Tableau 1-1. Classement de taille de cabine d’un avion selon la FAA [9].
≥71 Grand
30-70 Moyen
5-29 Petit
pas adaptée pour dégivrer les ailerons et gouvernes arrières ; les contraintes de
volume et de poids imposées en aéronautique rendent cette solution inadaptée pour
être utilisée au niveau des surfaces sensibles au givrage autres que les ailes de
l’avion,
interférence avec l’aérodynamique ; au moment de son activation, le pilote peut
sentir une perturbation dans le contrôle de la direction de l’aéronef,
problèmes de blocage mécanique ; il est arrivé que sous conditions de givrage
extrêmes les boudins restent bloqués à cause de la glace, soit en position initiale,
soit en position d’activation ; dans ces situations il est impossible de récupérer le
fonctionnement normal ; situation rare mais très dangereuse,
non viable pour les avions à turboréacteurs à double flux (turbine) ; quand on parle
des avions commerciaux, d’une taille au-delà de 100 tonnes, le système n’est plus
envisageable car le poids et la taille du compresseur et des connexions nécessaires
14
Chapitre I : Le dégivrage en aéronautique et ses technologies
Ces désavantages limitent l’utilisation des boudins pneumatiques à des avions privés de
propulsion à turbo moteur et moteur à pistons. Les tendances en termes d’avions de ligne
sont de concevoir des avions de tailles modérées à grandes [10], comme mis en évidence
sur la Figure 1-4. D’où la nécessité de trouver des solutions nouvelles.
Cette technologie démontre néanmoins qu’une action mécanique est très efficace pour
réaliser le travail de dégivrage. Plusieurs études sur la glace portent sur les effets des efforts
mécaniques sur les surfaces givrées [12]. Ces faits appuient la recherche de nouvelles
solutions hybrides avec action mécanique.
15
Chapitre I : Le dégivrage en aéronautique et ses technologies
C’est une solution très fiable même sous conditions givrantes sévères, l’idée a eu
ses origines dans les années 1920 [13], depuis l’apparition des avions à turboréacteurs à
double flux, et a été largement répandue vers les années 1940 [14]. On trouve cette
technologie dans pratiquement tous les avions commerciaux contemporains. La Figure 1-5
représente une vue transversale de l’aile dans laquelle la direction du flux d’air est montrée.
Figure 1-5. Système de circulation d’air chaud dans le bord d’attaque des ailes implanté sur des avions
Airbus, Dassault et Boeing [15] [16].
L’air chaud est conduit vers les tuyaux, surfaces et systèmes givrés. Le matériau du blindage
du bord d’attaque de l’aile est souvent un alliage d’aluminium ; les conduits d’air chaud en
titane ou acier inoxydable sont typiquement isolés avec de la fibre de verre, résistant au feu
et hautement isolant. En fonctionnement, la température externe des parties protégées par
le système antigivrage est de 60°C. La température maximale de sécurité est de 180°C du
coté intérieur et de 100°C du côté extérieur de l’aéronef [16]. La pression est aussi un
paramètre très surveillé dans ces systèmes.
16
Chapitre I : Le dégivrage en aéronautique et ses technologies
Bien que cette technique soit très performante pour le travail d’antigivrage, la solution reste :
hautement énergivore,
exclusivement adaptée aux avions à turbine.
Figure 1-6. Antigivrage par tapis électrothermique, pour le B787 de Boeing [15].
Les ETIPS peuvent intégrer une des deux technologies suivantes pour augmenter la
température : résistive ou par induction. La technologie résistive utilise une configuration
d’éléments résistifs qui vont chauffer via la circulation d’un courant continu, contrôlé
électroniquement et transférer la chaleur par conduction vers la surface à dégivrer. La
technologie inductive repose sur la circulation de courants induits pour chauffer directement
la surface à dégivrer. Le contrôle des champs électromagnétiques est également effectué
électroniquement. L’utilisation de cette technologie nécessite une modélisation rigoureuse,
les courant électriques et les charges électromagnétiques pouvant interférer avec les autres
dispositifs.
17
Chapitre I : Le dégivrage en aéronautique et ses technologies
ajout important de poids, dès lors que ce type de système doit couvrir l’intégralité de
la surface des ailes et des ailerons,
énergivores, ce sont des systèmes indépendants, dédiés, non intégrés aux autres
fonctions de l’avion qui absorbent une grande quantité d’énergie pour réaliser le
travail. Leur consommation énergétique n’est pas optimale,
besoin d’une source d’alimentation dédiée.
Les avions turbopropulsés, comme celui de la Figure 1-7, sont un exemple d’avion utilisant
une combinaison de systèmes électrothermiques et boudins pneumatiques pour la
protection contre le givre. De cette façon, on trouve un bon compromis entre la
consommation énergétique et l’efficacité de la protection antigivrage [6].
Figure 1-7. Avion de petite taille à usage particulier turbo propulsé avec protection antigivrage hybride.
18
Chapitre I : Le dégivrage en aéronautique et ses technologies
Figure 1-8. Configuration du dégivreur électromécanique conçu par Cox & Company, Inc [19].
Comme deuxième exemple, nous pouvons citer le brevet déposé en 2003 par la Goodrich
Corporation, compagnie spécialisée dans les systèmes de dégivrage et antigivrage. Cette
fois on peut se rendre compte, sur la Figure 1-9, qu’il s’agit d’une structure de dégivreur plus
compacte, qui est aussi fixée sous le blindage de l’aile. Il est prévu d’espacer les actionneurs
pour couvrir une surface plus importante de dégivrage.
19
Chapitre I : Le dégivrage en aéronautique et ses technologies
Figure 1-9. Schématiques du dégivreur électro expulsif attaché, conçu par Goodrich Corporation [20].
Le principe de ces technologies est d’injecter des courants de forte intensité et de sens
contraire, à travers les spires conductrices, de sorte à générer une force de répulsion entre
les spires, qui va venir s’appliquer à l’aile pour réaliser le dégivrage.
Au cours des dernières années, l’électronique n’a cessé d’évoluer et aujourd’hui les
systèmes électriques sont de taille plus réduite. Si en plus on cherche à repenser et
optimiser le mode de fonctionnement de ce type de solutions, il peut être intéressant de
reprendre cette technologie pour arriver à la conception d’un dispositif performant et
compatible avec l’avion plus électrique.
20
Chapitre I : Le dégivrage en aéronautique et ses technologies
De nombreuses solutions sont imaginées, mais elles nécessitent une validation pour
évaluer leur faisabilité, performance, consommation énergétique et intégrabilité. Les efforts
pour arriver à la conception de l’avion électrique poussent ces recherches vers l’avant. Dans
cette section, on présente une liste non exhaustive des technologies prometteuses qui sont
en voie de test et développement, parallèlement à la technologie étudiée dans cette thèse
(Figure 1-10).
21
Chapitre I : Le dégivrage en aéronautique et ses technologies
Peintures chauffantes. - L’idée est d’envelopper une série d’électrodes fines entre deux
couches de peinture de l’avion. Ces électrodes sont ensuite alimentées électriquement
conduisant à un échauffement par effet Joule [27]. Cela présenterait l’avantage d’un
système intégré dans la structure d’origine, avec un poids réduit.
Les recherches pointent donc vers l’intégration des technologies pour concevoir des
systèmes hybrides capables de fournir une somme d’avantages en essayant de minimiser
les inconvénients. Pour récapituler, le Tableau 1-2 permet de visualiser qualitativement les
22
Chapitre I : Le dégivrage en aéronautique et ses technologies
points forts et points faibles des principales technologies contre le givre actuellement
employées.
L’objectif de la thèse est donc d’expliciter les phénomènes en jeu, pour déterminer
les fonctions de sensibilité permettant d’ajuster la conception de la technologie. Le dispositif
de dégivrage imaginé par Zodiac, repose sur une action électromécanique. On posera les
bases de sa modélisation en vue de dégager des pistes vers l’optimisation de cette
technologie. Les détails de ce principe et le dispositif imaginé sont présentés dans le
paragraphe suivant.
23
Chapitre I : Le dégivrage en aéronautique et ses technologies
Pour respecter ces conditions, le principe choisi consiste à générer un effort mécanique à
partir d’un phénomène électromagnétique. Cet effort sera appliqué à la surface à dégivrer
pour induire la rupture de la pellicule de glace et l’éliminer complètement à l’aide de la
pression aérodynamique qui viendra naturellement la balayer. La Figure 1-11 montre la
procédure envisagée.
Figure 1-11. Principe d’expulsion de la glace proposé par l’équipe R&D de Zodiac Aerospace [31].
Le système est placé dans l’espace compris entre la structure rigide de l’aile de l’avion et le
blindage extérieur en aluminium qui la recouvre, au niveau des bords d’attaque. Autrement
dit, la structure de l’aile sera le support sur lequel le système va reposer et le blindage sera
la surface sur laquelle le travail mécanique sera appliqué. Il faut toutefois veiller à ne pas
affecter les caractéristiques mécaniques du blindage pour ne pas générer une déformation
permanente qui pourrait nuire à l’aérodynamique des ailes de l’avion.
24
Chapitre I : Le dégivrage en aéronautique et ses technologies
autour des conducteurs. Les conducteurs sont alors le siège d’efforts mécaniques, les forces
de Laplace, liées à l’interaction des courants électriques et du champ magnétique [32]. Les
efforts agissent de manière opposée sur chaque conducteur, selon la configuration illustrée
dans la Figure 1-12.
Figure 1-12. Disposition des pistes conductrices pour la génération des efforts mécaniques [33].
25
Chapitre I : Le dégivrage en aéronautique et ses technologies
intermittent maximum
-5 (JAR-25 Appendix C)
continuous maximum
(JAR-25 Appendix C)
-10
-15
-20
-25
-30
-35
-40
0 5000 10000 15000 20000 25000 30000 35000
PRESSURE ALTITUDE (ft)
Figure 1-13. Domaine des conditions givrantes rencontrées par un avion commercial selon FAR25 [31].
Figure 1-14. Taille des gouttes en fonction de la quantité d’eau liquide pour quatre températures [31].
26
Chapitre I : Le dégivrage en aéronautique et ses technologies
- verglas ou glaze ice : il s’agit d’un type de glace qui se forme quand la température
au moment de l’impact de l’eau sur l’aile est positive ; dans ce cas, une fraction d’eau
gèle lors de l’impact, le reste de la goutte ruisselle et gèle au long de sa trajectoire
formant un verre de glace dense et transparent,
- givre blanc ou rime ice : ce type de glace se forme quand la température ambiante
est négative au moment de l’impact des particules de glace contre le profil de l’aile ;
l’eau gèle alors instantanément et forme un givre blanc et opaque, peu dense.
Le verglas, quand la couche est peu épaisse, est plus facile à briser avec une action
mécanique, car il présente une rigidité plus prononcée. Le givre blanc a plus d’élasticité. Il
absorbe donc les efforts mécaniques ce qui le rend plus difficile à expulser par l’application
de forces [17].
Un démonstrateur a été conçu et testé par l’équipe Zodiac Aerospace en 2013. Les
informations présentées dans les paragraphes 1.3.c et 1.3.d ont été fournis et sont propriété
de la société Zodiac Aerospace.
Le dispositif construit par l’entreprise a été imaginé avec la forme d’un tapis
parallélépipédique d’épaisseur faible pour faciliter l’intégration. La Figure 1-15, représente
l’assemblage final du tapis expulse ayant servi pour les premiers essais en soufflerie de la
technologie. On y trouve la description des différentes couches constitutives. Comme
2Cette altitude est convenable pour diminuer la résistance de l’air, minimiser la consommation du carburant et
éviter les zones de haute précipitation.
27
Chapitre I : Le dégivrage en aéronautique et ses technologies
mentionné dans la description du principe physique, l’axe Z sera l’axe de référence pour
définir la direction de déplacement principale pour le dégivrage [31].
Figure 1-15. Assemblage du tapis expérimental de dégivrage développé par Zodiac [33].
Les pistes en rouge sont les conducteurs, disposés en forme de serpentin, qui vont se
déformer sous l’action des forces de Laplace. Les plots en polyamide auront le rôle de
ressorts mécaniques pour assurer le retour des éléments conducteurs à leur position initiale
après chaque activation. Le reste des éléments sont des matériaux de structure et de
conditionnement qui vont néanmoins avoir un impact sur le comportement mécanique du
système complet.
La source d’alimentation et les connectiques ne sont pas illustrées sur la Figure 1-15,
seulement l’actionneur. Le Tableau 1-3 contient la description plus détaillée et les
28
Chapitre I : Le dégivrage en aéronautique et ses technologies
Une lamelle latérale est l’élément unitaire d’une piste latérale et est
découpée dans une feuille de cuivre de 100µm. Elles permettent
24mm x 5mm x
C2 ou C4 Lamelle latérale l’ajustement de la valeur de la résistance finale du réseau. Elles ne Cuivre
100µm
participent que très peu au déplacement et sont donc considérées
comme immobile (pertes).
Constitué de deux plis d’isolants électriques eux mêmes percés de 30
trous de diamètre 2mm permettant le collage des plots de rappel sur
Polyimide 227mm x 84mm x
C3 (2) Ecran les bandes inter-piste. Il a pour fonction le rappel mécanique de la
(Kapton) 25µm
couche supérieur C5 et l’isolation électrique des deux réseaux
inférieur et supérieur.
(3) Support Feuille de polyimide découpée à la forme d’un parallélogramme avec
Polyimide 231mm x 89mm x
C1 ou C5 inférieur et un angle de 31°. Joue le rôle de support mécanique et d’isolant
(Kapton) 25µm
supérieur électrique vis-à-vis de l’extérieur.
L’inter-piste est l’espace compris entre deux pistes successives.
Dans le cas présent il est constant. Par extension, la bande inter-piste
(4) Bandes inter- Gomme
C1 ou C5 est le matériau résilient collé dans cette espace. Il garantit à la fois 250µm
piste élastomère
l’isolement électrique entre pistes et délimite leur jeu mécanique
latéral. Ce sont sur ces bandes que sont collés les plots de rappel.
Disposés en matrice de 3x10 et réalisés dans un morceau de tapis
élastomère sous forme de plot de 3mm de diamètre. Essentiellement
Gomme
C3 (5) Plots de rappel destinés à maintenir le réseau supérieur en position de repos. Ils 2 x 250µm
élastomère
agissent comme un ressort afin de ramener la partie supérieure de
l’actionneur dans son état initial.
Le motif ou topologie impose les emplacements de liaison des pistes
Motif / Topologie (connexion électrique), l’interpiste ainsi que la largeur de celles-ci. Il
détermine la longueur d’onde de la déformation mécanique.
29
Chapitre I : Le dégivrage en aéronautique et ses technologies
Le dispositif est vulcanisé dans l’intérêt de rapprocher les lignes de conducteurs pour
minimiser l’écart entre eux et maximiser l’amplitude des forces de Laplace générées, mais
aussi pour garantir sa robustesse face aux déformations consécutives. En conséquence,
les plots élastomères sont également comprimés. Sur la Figure 1-17 on peut apprécier le
sens d’application de la pression pendant la vulcanisation.
30
Chapitre I : Le dégivrage en aéronautique et ses technologies
La compression est réalisée sans vide d’air afin de ne pas réduire la déformation du
dispositif. La face en contact avec la face interne du blindage en aluminium correspond au
côté plat de la tôle de vulcanisation.
Les valeurs des épaisseurs après vulcanisation sont estimées, mais aucune mesure
n’a été réalisée après vulcanisation. Il faut aussi noter que ces informations sont valables
sous la condition qu’aucun élément ne s’est déplacé pendant la procédure de vulcanisation.
Tous les éléments élastomères sont collés avec de la colle néoprène. Un tissu
polyamide calandré recouvre les lamelles latérales pour assurer une meilleure isolation
électrique. La contribution de ces lamelles à la déformation est minime, elles ont
principalement pour but d’assurer la continuité du circuit électrique entre les lignes de
lamelles longitudinales.
Pour finaliser, le Tableau 1-4 résume les propriétés des matériaux utilisés pour la
fabrication de ce tapis démonstrateur.
Le tapis électro-expulse est installé sous une plaque en aluminium type 2024 de 0.5
mm d’épaisseur qui simule le blindage de l’avion [37]. Pendant les tests, le système n’est
pas chargé avec du givre. Il opère alors à vide. Les mesures du déplacement sont réalisées
avec un capteur laser aptes à mesurer le déplacement sur une position à la fois. Pour
31
Chapitre I : Le dégivrage en aéronautique et ses technologies
réaliser ces mesures, des points de référence ont été définis. Les coordonnées de ces points
A, B, C, D, E, qui se trouvent au milieu de la longueur du tapis, sont définis sur la Figure
1-18.
Figure 1-18. Plaque en aluminium et points de mesure, le tapis expulse est situé en-dessous [38].
Les points B et D sont alignés avec la position des lamelles conductrices du serpentin du
tapis. Les points A, C et E sont alignés avec la position des plots de rappel décrits sur la
Figure 1-15.
La Figure 1-19 contient les formes d’onde des courbes de tension et courant pour
un cycle d’activation d’une période estimée de T = 1,6 ms.
32
Chapitre I : Le dégivrage en aéronautique et ses technologies
Figure 1-19. Courant et tension d’alimentation pour une activation du tapis dégivrant [38].
La Figure 1-20 montre les résultats de mesure du déplacement au point B (0,-1), position
d’une ligne de plots de rappel, consécutif à un cycle d’alimentation. La Figure 1-21 montre
les résultats de mesure pour le déplacement au niveau du point central C (0,0)
correspondant à la position d’une lamelle du serpentin.
Figure 1-20. Déplacement en fonction du temps mesuré sur le point B pour un cycle d’activation [38].
33
Chapitre I : Le dégivrage en aéronautique et ses technologies
Figure 1-21. Déplacement en fonction du temps mesuré au point C pour un cycle d’activation [38].
34
Chapitre I : Le dégivrage en aéronautique et ses technologies
Ces tests ont permis d’illustrer le principe physique et d’obtenir des premiers
résultats pour démontrer la faisabilité de la technologie. Mais pour valider sa faisabilité et
35
Chapitre I : Le dégivrage en aéronautique et ses technologies
Celle-ci doit s’appuyer sur une description multiphysique afin d’identifier les
paramètres sensibles de conception adaptés aux travaux d’optimisation pour déterminer le
futur de la technologie.
A partir de ces constats s’impose le cahier des charges suivant (Tableau 1-5):
En plus de s’adapter au cahier des charges, l’actionneur devra respecter des critères
de masse, de volume et de durée de vie, être modulaire pour ne pas impacter le Mean Time
Between Failure ou taux de panne, de telle sorte que la solution reste intégrable.
36
Chapitre I : Le dégivrage en aéronautique et ses technologies
Dans cette thèse nous allons donc nous focaliser sur la modélisation d’un dispositif
de dégivrage et l’identification des paramètres de conception sensibles qui permettront
d’atteindre les objectifs imposés par le cahier des charges.
I.5 Conclusion
La solution a été décrite dans la deuxième partie du présent chapitre, elle a été
évaluée au préalable par le partenaire industrielle Zodiac à travers la construction et le test
d’un démonstrateur. Les données expérimentales issues de ce prototype appuient les
travaux d’analyse et de modélisation dans une perspective d’optimisation du dispositif. On
présume qu’avec une correcte conception de cette solution, la quantité d’énergie utilisée
peut être optimisée.
Cette démarche est intéressante dans la mesure où elle devrait permettre d’assurer
l’obtention d’un modèle le plus simple possible, avec le minimum de paramètres nécessaires
pour représenter les phénomènes et la construction du dispositif, sans pour autant être
37
Chapitre I : Le dégivrage en aéronautique et ses technologies
38
Chapitre I : Le dégivrage en aéronautique et ses technologies
39
Chapitre I : Le dégivrage en aéronautique et ses technologies
[21] Anne de Becco, “Rapport annuel Zodiac Aerospace France 2012-2013,” Plaisir France,
Communication, 2013.
[22] P. Kim, T.-S. Wong, and & al., “Liquid-infused Nanostructured Surfaces with Extreme
Anti-ice and Anti-Frost Performance,” presented at the ACS Nano, 2012, pp. 6569–
6577.
[23] Vahid Hejazi, Konstantin Sobolev, and Michael Nosonovsky, “From superhydrophobicity
to icephobicity: forces and interaction analysis,” Nat. Res. J., vol. Scientific Reports 3,
Jul. 2013.
[24] [Link], [Link], [Link], [Link], and [Link], “Understanding the effect
of superhydrophobic coatings on energy reduction in anti-icing systems,” Cold Regions
Science and Technology, vol. 67, pp. 59–67, Feb-2011.
[25] Yoshinori Furukawa and John S. Wettlaufer, “Snow and ice crystals,” Physics Today,
vol. 60, no. 12, p. 2, 2007.
[26] S.V. Venna and Y.-J. Lin, “Mechatronic Development of Self-Actuating In-Flight Deicing
Structures,” in IEEE/ASME Transactions on Mechatronics, 2006, vol. 11, pp. 585–592.
[27] J. Henrion, M. Claire, “Peinture chauffante amliorée, Brevet A1, Published on June
2014, 19th.,” WO2014091161 A1, 19-Jun-2014.
[28] Strehlow, R. and Moser, R., “Capitalizing on the Increased Flexibility that Comes from
High Power Density Electrothermal Deicing,” Aerospace Technology Conference and
Exposition, 2009.
[29] T. Wright, “Electro-mechanical Deicing,” Air & Space Magazine, Mar-2004. .
[30] Haslim, “IN/ARZIN/ARZ: Réévaluation d’un brevet NASA,” 44690 69035335.
[31] S. L. Garrec, E. Belot, and A. Delehelle, “Spécification d’un Actionneur
Electromécanique,” Zodiac Aerospace, Intertechnique 7N–17079, Oct. 2013.
[32] J. D. Jackson, Classical Electrodynamics. New York, USA: John Wiley & Sons, 1962.
[33] E. Belot, “Dossier technique de définition d’un actionneur expulse,” Zodiac Aerospace,
Plaisir, France, 7N–17211, Apr. 2013.
[34] Federal Aviation Administration, “eCFR Regulations.” [Online]. Available:
[Link]
idx?SID=a8bb87823ac0ab79d7adae35710d7c44&mc=true&node=pt14.1.25&rgn=div5
.
[35] C. Lange, J. Baldzuhn, S. Fink, R. Heller, M. Hollik, and W. H. Fietz, “Paschen Problems
in Large Coil Systems,” presented at the IEEE Transactions on applied
superconductivity, 2012, vol. 22.
[36] G. R. Govinda Raju and R. Hackam, “Note on Paschen Law and the Similarity Theorem
at the minimum breakdown voltage,” IEEE transactions on Plasma Science, pp. 63–66,
1974.
[37] ALCOA, “ALLOY 2024 SHEET AND PLATE.” Alcoa Mill Products Inc.
[38] Eric Belot, “Réponse_BA_HY_NRC2012_EI1_@[Link].” Zodiac Services Europe -
Groupe Zodiac, 04-Nov-2013.
[39] Zodiac Aerospace France, BA_HY_EI1_LargeAngle Animation du résultat des mesures
de déplacement pour le demonstrateur du actionneur expulse. 2013.
40
Chapitre 2 :
Analyse du dégivreur et modèle
électromagnétique
Chapitre 2 : Analyse du dégivreur et modèle électromagnétique
II.3.b Modèle électromagnétique retenu et outil de calcul retenu pour l’approche quasi
statique 54
II.1 Introduction
Deux approches ont été abordées. La première est une approche harmonique pour
prendre en compte la dynamique du signal électrique d’alimentation, qui va être
abandonnée pour des raisons qui seront expliquées. Une deuxième approche, de type
calcul en statique, a été retenue. Celle-ci a été intégrée dans une plateforme adaptée pour
un calcul quasi-statique qui permettra la prise en compte partielle de la dynamique des
phénomènes électromagnétiques. Cette plateforme acceptera, dans un deuxième temps,
l’intégration du modèle multiphysique complet final.
43
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-2. Vue latérale, plan YZ ; du serpentin en cuivre à modéliser et sa vue en détail.
Ces serpentins sont alimentés par une impulsion de courant d’amplitude suffisamment
importante pour que la force de répulsion engendrée entre les conducteurs puisse les
déformer, et en conséquence, déformer aussi la surface à dégivrer. La déformation
mécanique doit rester dans le domaine élastique des matériaux, de façon à ce que la glace
accrétée sur le bord d’attaque de l’aile de l’avion puisse être expulsée à chaque activation
du dégivreur [1].
44
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-3. Vue en section d’une ligne conductrice avec les lamelles fusionnées.
45
Chapitre II : Analyse du dégivreur et modèle électromagnétique
La Figure 2-5 indique la direction du courant (signalée avec un point ou une croix pour des
sens contraires de circulation du courant) ainsi qu’un exemple de la force et de l’induction
magnétique produite par l’interaction entre deux conducteurs quelconques dans le plan YZ.
Figure 2-5. Coupe transversale du serpentin sur un plan YZ (pour une valeur d’abscisse constante).
Selon la convention d’axes retenue, il ressort que c’est la composante Y du champ B qui
va créer la composante Z de la force qui induira la déformation de la structure (ANNEXE
A). L’induction magnétique totale vue pour un conducteur résulte de la contribution de tous
les courants associés aux conducteurs selon la disposition de la Figure 2-5.
𝜇0 ∗ 𝐼
𝐵= (2-1)
2𝜋 ∗ 𝑟
46
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-6. Modèle de conducteur filiforme infiniment long, à section nulle a) et à section circulaire
non nulle b).
𝜇𝑜 ∙ 𝐼 ∙ 𝑟
( 2
) 𝑠𝑖 0 ≤ 𝑟 ≤ 𝑅
𝐵(𝑟) = { 2𝜋 ∙ 𝑅 (2-2)
𝜇𝑜 ∙ 𝐼
( ) 𝑝𝑜𝑢𝑟 𝑟 > 𝑅
2𝜋 ∙ 𝑟
Cependant, la section des conducteurs du dispositif étudié n’est pas circulaire mais
rectangulaire, et la densité de courant n’est pas uniforme du fait des effets peau et de
proximité. L’approche analytique s’avère donc insuffisante. Pour la suite des modélisations
électromagnétiques un modèle semi-analytique est mise en œuvre. La méthode de
résolution PEEC (Partial Element Equivalent Circuit), ou résolution par intégration
d’éléments partiels d’un circuit électrique, est utilisée. Cette modélisation repose sur une
résolution harmonique pour le calcul rapide de circuits impédants [2].
47
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-7. Induction magnétique pour un conducteur filiforme et à section rectangulaire et les
simulations InCa3D.
Par ailleurs, la simulation InCa3D présente les résultats seulement pour le premier
harmonique (f= 625 Hz) du courant d’alimentation. La prise en compte de l’effet peau et de
proximité a été effectuée pour cette fréquence uniquement. La comparaison des ordres de
grandeur entre les deux approches n’a pas d’intérêt, cependant son évolution permettrait
de valider la direction et le sens du champ magnétique de ces deux approches.
48
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Dans ce paragraphe, la force de Laplace entre les deux serpentins conducteurs est
évaluée. La force électromécanique permet le travail mécanique à l’origine du dégivrage.
Elle représente l’effort fourni par le système, et sera une donnée d’entrée pour le
développement du modèle mécanique.
49
Chapitre II : Analyse du dégivreur et modèle électromagnétique
𝑑𝐹⃗ = 𝐼 ∗ ⃗⃗⃗⃗ ⃗⃗
𝑑𝑙 ∧ 𝐵 (2-3)
50
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-9. Géométrie simulée sur InCa3D, vue isométrique et plan XZ.
Le Tableau 2-1 indique les données d’entrée de cette modélisation ainsi que le résultat
obtenu. Ce résultat correspond à la force moyenne par conducteur pour un circuit alimenté
avec un signal de courant DC.
Tableau 2-1. Valeurs d’entrée pour les paramètres géométriques et physiques des premières
simulations pour la force.
51
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Décomposition harmonique
Comme précédemment indiqué, InCa3D réalise une modélisation harmonique.
Aussi, pour travailler dans le domaine fréquentiel la transformée de Fourier rapide (FFT) du
signal du courant est effectuée. La procédure appliquée est la suivante : chaque
harmonique de courant est utilisé comme courant d’alimentation dans la modélisation
PEEC et le principe de superposition est utilisé pour reconstruire la réponse dans sa totalité.
Une transformée inverse (FFT-1) est ensuite appliquée pour revenir à une description
temporelle.
52
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Les composantes continues (en rouge) et pulsatoire (en bleu) sont bien calculées
par InCa3D, mais les autres termes (en vert), faisant intervenir deux harmoniques différents
simultanément, ne sont pas pris en compte dans InCa3D.
𝑑𝐹(𝑡) = 𝐼𝑜 2 + 𝐼1 2 + 𝐼2 2
+𝐼1 2 cos(2𝜔1 𝑡 + 2𝜑1 ) + 𝐼2 2 cos(2𝜔2 𝑡 + 2𝜑2 )
+2√2𝐼𝑜 𝐼1 cos(𝜔1 𝑡 + 𝜑1 ) + 2√2𝐼𝑜 𝐼2 cos(𝜔2 𝑡 + 𝜑2 )
+𝟐𝑰𝟏 𝑰𝟐 [𝐜𝐨𝐬((𝝎𝟏 + 𝝎𝟐 )𝒕 + 𝝋𝟏 + 𝝋𝟐 ) + 𝐜𝐨𝐬((𝝎𝟏 − 𝝎𝟐 )𝒕 + 𝝋𝟏 − 𝝋𝟐 )] (2-5)
3
𝜇0 ∙𝑑𝑙
Le coefficient , devant chacune des composantes, a été mis à 1 pour plus de simplicité et clarté dans
2𝜋∙𝑟
la comparaison.
53
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Pour des conducteurs alimentés par des courants de sens contraire et constitués
d’un fondamental et de deux harmoniques, il est possible d’apprécier le terme d’interaction
entre harmoniques. La Figure 2-10 indique clairement que cette interaction entre
harmoniques ne peut pas être négligée dans le calcul de la force magnétique. Le modèle
PEEC intégré dans l’outil InCa3D ne peut pas être retenu pour la modélisation
électromagnétique. Une solution pour contourner cette difficulté repose sur l’utilisation
d’une modélisation statique dans un schéma itératif dans le temps.
Dans un second temps, il sera nécessaire de trouver une solution pour reconsidérer
la dépendance au temps. Comme il est indiqué dans le chapitre 4, une méthode de
résolution itérative, pas à pas dans le temps, a été mise en œuvre.
54
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Une méthode de calcul numérique est proposée pour la résolution des équations
mathématiques pour le système. La résolution du modèle sera en mode quasi-statique.
55
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Le model retenu est créé à l’aide des logiciels de simulations dédiés énoncés
antérieurement et ensuite décrits plus en détail.
Mise en œuvre
Le modèle électromagnétique retenu est construit à partir de la définition
géométrique de deux conducteurs, motif élémentaire extrait du serpentin original (§ II.3
Hypothèse 3), et réalisée avec le logiciel McMMems.
56
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Pour cette première étude, la description des conducteurs reste très simple : le changement
de position généré par la déformation mécanique n’est pas pris en compte, un seul
conducteur massif non déformable est défini. Les calculs ont été comparés à ceux de
InCa3D dans le cas d’un courant d’amplitude constante et validés avec des pourcentages
d’erreur négligeables pour le calcul de la force (erreur<0.1%).
Scénario 1
Le Tableau 2-2 contient les données d’entrée du premier scénario. Les résultats
pour la force par conducteur en fonction de la variation du courant DC d’alimentation du
circuit sont présentés sur la Figure 2-12.
Tableau 2-2. Valeurs d’entrée pour les paramètres géométriques et physiques du 1 er scenario.
57
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-12. Force de Laplace produite sur les conducteurs en fonction du courant d'alimentation.
La convergence des résultats obtenus avec les deux logiciels, InCa3D et CADES Calculator
avec le composant *.icar valide la modélisation développée. Comme attendu, l’amplitude
de la force de Laplace est proportionnelle au produit carré de l’intensité du courant
d’alimentation, pendant que le sens de la force est contraire entre les deux conducteurs.
Une première piste pour l’optimisation de la conversion d’énergie serait alors de réduire
l’impédance de notre dispositif pour avoir un courant maximal.
Scénario 2
Le Tableau 2-3 indique les données d’entrée du deuxième scénario pour la
géométrie du cas étudié. La force par conducteur en fonction de l’écart entre conducteurs
a été calculée pour 3 valeurs du courant d’alimentation. Les résultats pour un courant de I
= 500 A, I = 3500 A et I = 7000 A, sont représentés sur les Figure 2-13, Figure 2-14 et
Figure 2-15 respectivement.
Tableau 2-3. Valeurs d’entrée pour les paramètres géométriques et physiques du 2 ème scenario.
58
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-13. Force de Laplace produite sur les conducteurs en fonction de l’écart (avec I = 500 A).
Figure 2-14. Force de Laplace produite sur les conducteurs en fonction de l’écart (avec I = 3500 A).
59
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-15. Force de Laplace produite sur les conducteurs en fonction de l’écart (avec I = 7000 A).
Scénario 3
Le Tableau 2-4 indique les données d’entrée du troisième scénario. Comme pour le
scenario 2, 3 valeurs du courant d’alimentation ont été testées. Les résultats concernant la
force par conducteur en fonction de la variation de l’épaisseur des conducteurs en cuivre
sont représentés sur les Figure 2-16, Figure 2-17 et Figure 2-18.
Tableau 2-4. Valeurs d’entrée pour les paramètres géométriques et physiques du 3ème scenario.
60
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-16. Force de Laplace produite sur les conducteurs en fonction de l’épaisseur (I = 500 A).
Figure 2-17. Force de Laplace produite sur les conducteurs en fonction de l’épaisseur (I = 3500 A).
61
Chapitre II : Analyse du dégivreur et modèle électromagnétique
Figure 2-18. Force de Laplace produite sur les conducteurs en fonction de l’épaisseur (I = 7000 A).
Conclusions
Le dispositif étudié est alimenté par le courant de décharge d’un condensateur. La
conséquence en est la création d’efforts de répulsion électromagnétiques entre les
différents conducteurs qui le constituent. Un modèle qui permet d’évaluer les efforts
électromagnétiques entre les conducteurs de ce dispositif pour une valeur de courant
donnée a été mis en place.
62
Chapitre II : Analyse du dégivreur et modèle électromagnétique
la géométrie fidèle du dispositif est reprise pour générer un modèle statique qui a été évalué
et a permis de valider les résultats obtenus avec les logiciels expérimentaux.
63
Chapitre II : Analyse du dégivreur et modèle électromagnétique
[1] E. Belot, “Dossier technique de définition d’un actionneur expulse,” Zodiac Aerospace,
Plaisir, France, 7N–17211, Apr. 2013.
[2] E. Clavel, “Vers un Outil de Conception de Câblage : Le Logiciel InCa.,” Electric power,
Institut National Polytechnique de Grenoble, France, 1996.
[3] J. D. Jackson, Classical Electrodynamics. New York, USA: John Wiley & Sons, 1962.
[4] H. L. RAKOTOARISON, “Méthode et outil de génération automatique de modèle pour
l’optimisation fortement contrainte des microsystèmes magnétiques,” Université Joseph
Fourier, Grenoble France, 2007.
[5] H. L. Rakotoarison, V. Ardon, O. Chadebec, B. Delinchant, and S. Guérin, “Formal
Sensitivity Computation of Magnetic Moment Method.,” IEEE Transac- tions on
Magnetics, Institute of Electrical and Electronics Engineers, pp. 1014–1017, Jun-2008.
[6] H. L. Rakotoarison, F. Wurtz, B. Stephanek, and O. Cugat, “Integration of the volume
segmentation method into an optimization process: application to the sizing of a micro-
actuator for deformable mirrors,” presented at the IEEE Transactions on Magnetics, vol.
42, pp. 1163–1166.
[7] P. PHAM-QUANG, “Modélisation magnéto-mécanique d’un nano commutateur.
Optimisation sous contraintes de fiabilité par dérivation automatique des programmes
en Java,” Ph.D. Génie Electrique, Université de Grenoble, Grenoble France, 2011.
64
Chapitre 3 :
Modélisation mécanique
Chapitre 3 : Modélisation mécanique
III.1 Introduction ........................................................................................................... 67
III.1 Introduction
Dans cette partie, le cadre de l’étude et les notions de mécanique des milieux
continus et de mécanique des structures sont rappelés. Ces notions sont la base pour
l’étude de la déformation de la plaque à modéliser. Sous les hypothèses de déformation
67
Chapitre III : Modélisation mécanique
domaine 𝛺.
Figure 3-1. Définition du domaine, frontières, efforts et déplacements d’un milieu continu.
La mécanique des structures est la mécanique des solides de dimensions finies. Elle
découle de la mécanique des milieux continus et exploite une caractéristique essentielle
des structures : une dimension au moins est faible devant les autres. En d’autres termes,
les trois dimensions d’une structure ne sont pas du même ordre de grandeur. Les
géométries les plus communes sont les poutres (1D), les plaques et les coques (2D). La
contrepartie à cette simplification est une perte de généralité : la mécanique des
structures donne une représentation satisfaisante sous conditions qu’un certain nombre
d’hypothèses soient satisfaites. La définition de ces hypothèses cinématiques est la
68
Chapitre III : Modélisation mécanique
principale difficulté. Elle doit être suffisamment riche pour décrire tous les phénomènes
observés et assez simple pour permettre des résolutions analytiques.
Pour toute étude, des milieux continus homogènes et isotropes sont considérés.
Le tenseur de Green-Lagrange définit les changements relatifs de longueur du domaine 𝛺
par rapport à un état de référence. Dans le cas de petites perturbations avec de petits
déplacements, ce tenseur définit une relation linaire entre les déplacements et les
déformations :
1 𝜕𝑢𝑖 𝜕𝑢𝑗
𝑆𝑖𝑗 = ( + ) 𝑖, 𝑗 = 1,2,3. (3-1)
2 𝜕𝑥𝑗 𝜕𝑥𝑖
𝜎 = 𝐶: 𝑆 , (3-2)
Voigt) sous la forme d’une matrice 6x6 où les directions représentent les directions de
déformation. Dans le cas des milieux isotropes, le tenseur d’élasticité est défini à partir
des coefficients de Lamé 𝜆∗ et 𝜇∗ :
𝐶𝑖𝑗𝑘𝑙 = 𝜆∗ 𝛿𝑖𝑗 𝛿𝑘𝑙 + 2𝜇 ∗ 𝛿𝑖𝑘 𝛿𝑗𝑙 , (3-3)
𝜐∙𝐸 𝐸
où : 𝜆∗ = 𝜇∗ =
(1+𝜐)(1−2𝜐) 2(1+𝜐)
69
Chapitre III : Modélisation mécanique
𝑅 = (∫ (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗𝑑𝑠)
𝑒1 𝑒1 + (∫ (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗𝑑𝑠)
⃗⃗⃗⃗ 𝑒2 𝑒2 + (∫ (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗𝑑𝑠)
⃗⃗⃗⃗ 𝑒3 𝑒3
⃗⃗⃗⃗ (3-4)
∑ ∑ ∑
N Vy Vz
(𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗
𝑒1
𝑜
𝑀 = ∫ (𝑦) ⋀ (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗
𝑒2 𝑑𝑠
∑ 𝑧
= (∫ (𝑦 (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗
𝑒3 − 𝑧 (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗)
𝑒2 𝑑𝑠) ⃗⃗⃗⃗
𝑒1 + (∫ (𝑧 (𝜎 ∙ 𝑛⃗⃗) ∙ ⃗⃗⃗⃗)
𝑒1 𝑑𝑠) ⃗⃗⃗⃗
𝑒2
∑ ∑
Mx My
Mz
70
Chapitre III : Modélisation mécanique
Figure 3-2. Exemple d’un tronçon de poutre et représentation des projections de la force de réaction.
𝑑 2 𝑢(𝑥⃗, 𝑡)
𝜌 = 𝜌𝑓⃗(𝑥⃗, 𝑡) + ⃗⃗⃗⃗⃗⃗⃗
𝑑𝑖𝑣 𝜎(𝑥⃗, 𝑡) ∀ 𝑥⃗ ∈ ℜ (3-6)
𝑑𝑡 2
Cette relation d’équilibre est complétée par des conditions aux limites cinématiques
(Equation 3-7) et des contraintes (Equation 3-8) :
⃗⃗(𝑥⃗, 𝑡) = ⃗⃗⃗⃗⃗
𝑢 𝑢 𝑠 (𝑥⃗, 𝑡) ∀ 𝑥⃗ ∈ 𝜕Ω𝑢 (3-7)
Figure 3-3. Motif élémentaire pour l’analyse de poutres, réduction du domaine d’étude par symétrie ¼.
71
Chapitre III : Modélisation mécanique
Une poutre est un élément de structure allongé, de section réduite, avec un centre
de gravité G sur la longueur, soumise à des conditions à ses extrémités. La section de la
⃗⃗ et l’état de
poutre est petite face à sa longueur, de sorte que le champ de déplacement 𝑈
72
Chapitre III : Modélisation mécanique
𝑑𝑣(𝑥)
𝑢𝑚 (𝑥, 𝑦) = 𝑢(𝑥) − 𝑦 (3-12)
𝑑𝑥
𝑑𝑣(𝑥)
avec : = 𝜃(𝑥)
𝑑𝑥
73
Chapitre III : Modélisation mécanique
Dans le cadre de cette cinématique, les relations d’équilibre (Equation 3-6 à 3-8) que
doivent vérifier le déplacement de tout point 𝑚 d’une section, se réduisent à l’équation
d’Euler-Bernoulli4 [2].
𝑑2 𝑑2 𝑣(𝑥)
(𝐼 ∙ 𝐸 ∙ ) = 𝑝𝑦 (𝑥) (3-14)
𝑑𝑥 2 𝑧 𝑑𝑥 2
𝐼𝑧 = ∫ 𝑦 2 𝑑𝑠 ( 3-15)
𝑠
𝑑4 𝑣(𝑥)
𝐼𝑧 ∙ 𝐸 ∙ = 𝑝𝑦 (𝑥) (3-16)
𝑑𝑥 4
𝑑 𝑑 2 𝑣(𝑥)
(𝐼𝑧 ∙ 𝐸 )| = 𝑇 (3-17)
𝑑𝑥 𝑑𝑥 2 𝑥 0
𝑑2 𝑣(𝑥)
𝐼𝑧 ∙ 𝐸 | =𝑀 (3-18)
𝑑𝑥 2 𝑥
0
𝑑𝑣(𝑥)
| = 𝜃(𝑥0 ) (3-19)
𝑑𝑥 𝑥0
avec 𝑥0 = {0, 2}
4 L’étude de cette équation d’équilibre repose sur le fait que la contribution du terme de cisaillement est
inférieur au terme de rotation 𝜃(𝑥) lors du calcul de la flèche 𝑣(𝑥) et peut donc être négligé.
74
Chapitre III : Modélisation mécanique
⃗⃗1 = 𝐹𝑑1
𝑢 = 𝑢𝑑 𝑜u 𝑛⃗⃗ ∙ 𝑁 (3-20)
⃗⃗2 = 𝐹𝑑2
𝑣 = 𝑣𝑑 𝑜u 𝑛⃗⃗ ∙ 𝑁 (3-21)
⃗⃗⃗ = 𝑀𝑑𝑧
𝑣1𝑥 = 𝜃𝑑 𝑜u 𝑛⃗⃗ ∙ 𝑀 (3-23)
Le Tableau 3-1 définit les conditions aux limites associées à ces deux types de
comportement.
Tableau 3-1. Traduction des conditions aux limites au niveau des grandeurs mécaniques.
Le traitement de l’équation 3-13 tenant compte de ces conditions aux limites respectives
pour ces deux cas, conduit aux solutions analytiques du déplacement pour une poutre bi-
appuyée (Equation 3-24) et bi-encastrée (Equation 3-25) :
𝑃𝑦
v(𝑥) = (𝑥 4 − 2𝐿𝑥 3 + 𝐿3 𝑥) (3-24)
24 ∙ 𝐸 ∗ 𝐼𝑧
1
v(𝑥) = (𝑥 4 − 𝐿𝑥 3 + 𝐿2 𝑥 2 ) (3-25)
24 ∙ 𝐸 ∗ 𝐼𝑧
75
Chapitre III : Modélisation mécanique
La Figure 3-5 illustre le déplacement d’une poutre de longueur unitaire pour ces
deux types de conditions aux limites.
Figure 3-5. Déplacement analytique pour un chargement uniforme d’une poutre bi-appuyée ou bi-
encastrée.
76
Chapitre III : Modélisation mécanique
Figure 3-6. Représentation graphique des modèles d’une poutre discrétisée en n-éléments, appui et
encastrement.
a) b)
c) d)
Figure 3-7. Déformation pour une poutre bi-encastrée et une poutre bi-appuyée, pour un seul élément
poutre a) et avec une discrétisation en 3 éléments b), avec le modèle discrétisé pour un chargement
non uniforme et des conditions aux limites d’appui c) ou d’encastrement d).
77
Chapitre III : Modélisation mécanique
Sur la Figure 3-7 c) et d) on observe les courbes de déformation obtenues avec des
chargements non-uniformes pour les deux types de condition aux limites considérés, on
constate une continuité des courbes liée à une correcte prise en compte des conditions
de continuité.
Figure 3-8. Section de la poutre 3D, segment de blindage en aluminium (suffixe alu, i = 1) et
conducteur longitudinal en cuivre (suffixe cui, i=2).
𝐸𝑖
𝑌𝑒𝑞 𝑖 = 𝑌𝑖 (3-26)
𝐸𝑟𝑒𝑓
78
Chapitre III : Modélisation mécanique
∑𝑛𝑖=1 𝑍𝑝 𝑖 ∙ 𝑌𝑒𝑞 𝑖 ∙ 𝑍𝑖
𝑍𝑛 = (3-27)
∑𝑛𝑖=1 𝑌𝑒𝑞 𝑖 ∙ 𝑍𝑖
𝐼𝑛 𝑖 = 𝐼𝑧 𝑖 + 𝑆𝑖 ∙ 𝑑𝑖2 (3-28)
avec :
𝑌𝑒𝑞 𝑖 ∗ 𝑍𝑖3
𝐼𝑧 𝑖 = (3-29)
12
𝑑𝑖 = |𝑍𝑛 − 𝑍𝑝 𝑖 | (3-30)
𝑆𝑖 = 𝑌𝑒𝑞 𝑖 ∗ 𝑍𝑖 (3-31)
79
Chapitre III : Modélisation mécanique
Figure 3-9. Décomposition des conducteurs du modèle en deux lamelles par ligne.
Le découpage est voulu par le partenaire Zodiac pour les raisons suivantes :
Dans le but d’identifier si le contact a lieu, trois cas sont considérés par la
suite (Tableau 3-3).
Tableau 3-3. Description des trois cas comparatifs pour l’étude du feuilletage des conducteurs.
80
Chapitre III : Modélisation mécanique
Pour l’évolution de ces trois cas, les données géométriques sont indiquées dans le
Tableau 3-4. Les données des matériaux sont celles du Tableau 3-2.
Tableau 3-4. Données géométriques pour l’analyse de l’impact du feuilletage des conducteurs.
Données
𝐼 courant d’alimentation par couche/spire/niveau 5800 A
Epaisseur du conducteur massif 150 µm
Epaisseur d’une lamelle 75 µm
Ecart entre les deux conducteurs massifs ou entre les deux 125 µm
lamelles intérieures
Brasure (écart entre les deux lamelles du même conducteur) 25 µm
Longueur 223 mm
Largeur 3 mm
Epaisseur du blindage en aluminium 0,5 mm
81
Chapitre III : Modélisation mécanique
Le Tableau 3-5 recueille les résultats du bilan des efforts pour les trois cas,
calculés à l’aide du modèle électromagnétique final construit sur les logiciels McMMems
et CADES.
Tableau 3-5. Cas étudiés pour la décomposition des conducteurs en multiples lamelles.
82
Chapitre III : Modélisation mécanique
a)
83
Chapitre III : Modélisation mécanique
b)
Figure 3-10. Mise en évidence d’une interférence entre les conducteurs C2 et C3, déplacement avec
conditions aux limites du type appui a) et du type encastrement b).
Les résultats montrent que l’effort appliqué sur la lamelle intérieure, conducteur
C2, est suffisant pour qu’elle entre en contact avec la lamelle extérieure, conducteur C3
(Figure 3-10). Le déplacement calculé pour la lamelle supérieure interne est plus
important que celui de la lamelle supérieure externe. Il n’est physiquement pas possible
que le déplacement de la lamelle intérieure soit plus important que celui de la lamelle
extérieure. Une contrainte est manquante pour définir l’interaction mécanique entre les
lamelles des conducteurs à partir du moment où elles entrent en contact.
On observe aussi que les déplacements obtenus avec ce modèle sont de l’ordre
du mètre; très éloignés des déplacements mesurés lors des expériences menées par le
partenaire Zodiac donnant une valeur approchée du millimètre. Ces valeurs témoignent
des limitations de ce modèle et de la difficulté à prendre en compte les conditions aux
limites réellement en jeu dans le fonctionnement du dispositif.
Notons finalement qu’il n’y a pas d’isolation électrique (autre que l’air) entre les
deux lamelles. Il est difficile de connaitre les conséquences électriques et mécaniques de
ce contact car pour l’instant il n’est pas possible de déterminer le moment précis où le
contact s’établit.
La modélisation traitée ici repose sur la théorie des poutres qui favorise une
direction principale. Bien que les résultats ne soient pas exploitables, d’un point de vue
quantitatif liés aux hypothèses sous-jacentes à cette approche, cette première approche
simpliste permet une analyse qualitative du comportement du dispositif.
84
Chapitre III : Modélisation mécanique
Figure 3-11. Instantané du film fourni par ZODIAC montrant le déplacement de la structure à l’instant
t = 0,185 ms, vmax(x,y) = 0,35 mm [3] .
Ceci nous conduit à adopter pour la suite de la modélisation une approche de type
plaque et la présence d’efforts de rappel imposés par les plots élastomères.
Une plaque est un corps solide dont une des dimensions, l’épaisseur, est petite
devant les autres. Cela correspond, typiquement, à un rapport de la taille caractéristique
𝐿
de la plaque sur l’épaisseur ℎ
> 5. La surface de la plaque, ou feuillet moyen, est plane.
Comme dans le cas des poutres, les grandeurs ne sont plus définies en 3D, mais selon
l’épaisseur et le plan de la plaque. On définit alors des caractéristiques des forces par
unité de longueur (des contraintes intégrées sur l’épaisseur de la plaque).
85
Chapitre III : Modélisation mécanique
Par la suite, on considère que l’épaisseur de la plaque est supposée très petite
devant la taille caractéristique de la plaque. Il en découle que le cisaillement transverse
est négligé et les hypothèses suivantes sont ajoutées :
- le feuillet moyen ne subit pas de déformation dans son plan, seul le déplacement
transversal noté 𝑤(𝑥, 𝑦) des points du feuillet moyen est considéré,
- les sections normales au feuillet moyen restent normales lors de la déformation
(du fait que le cisaillement est négligé),
- hypothèse des petites déformations.
𝑢𝑚 𝑢(𝑥, 𝑦) + 𝑧 𝜃𝑦 (𝑥, 𝑦)
𝑢 𝑣
⃗⃗𝑚 = ( 𝑚 ) = (𝑣(𝑥, 𝑦) − 𝑧 𝜃𝑥 (𝑥, 𝑦)) (3-33)
𝑤𝑚 𝑤(𝑥, 𝑦)
𝜕𝑤 𝜕𝑤
avec : 𝜃𝑥 = et 𝜃𝑦 = − 𝜕𝑥
𝜕𝑦
membrane :
𝜕𝑁𝛼𝛽 (𝑥𝑚 )
+ 𝑝𝛼 (𝑥𝑚 ) = 0 (3-34)
𝜕𝛽
flexion :
𝜕 2 𝑀𝛼𝛽 (𝑥𝑚 )
+ 𝑝(𝑥𝑚 ) = 0 (3-35)
𝜕𝛼 𝜕𝛽
avec (𝛼, 𝛽) = {𝑥, 𝑦}, 𝑝𝛼 (𝑥𝑚 ) les efforts répartis dans le plan de la plaque, 𝑝(𝑥𝑚 ) la
pression et les relations de comportement pour une plaque orthotrope linéaire :
86
Chapitre III : Modélisation mécanique
𝐸𝑥 ∙ ℎ
∙ (𝜀𝑥𝑥 + 𝜐𝑦𝑥 + 𝜀𝑦𝑦 )
ℎ/2 1 − 𝜐𝑥𝑦 𝜐𝑦𝑥
𝑁𝛼𝛽 = ∫ 𝜎𝛼𝛽 𝑑𝑧 = 𝐸𝑦 ∙ ℎ (3-36)
−ℎ/2 ∙ (𝜀𝑦𝑦 + 𝜐𝑥𝑦 + 𝜀𝑥𝑥 )
1 − 𝜐𝑥𝑦 𝜐𝑦𝑥
{ 2 ∙ 𝐺𝑥𝑦 ∙ ℎ ∙ 𝜀𝑥𝑦
𝐸𝑥 ∙ ℎ3
∙ (𝑘𝑥𝑥 + 𝜐𝑦𝑥 + 𝑘𝑦𝑦 )
12(1 − 𝜐𝑥𝑦 𝜐𝑦𝑥 )
ℎ/2
𝐸𝑦 ∙ ℎ3
𝑀𝛼𝛽 = ∫ 𝑧 ∙ 𝜎𝛼𝛽 𝑑𝑧 = ∙ (𝑘𝑦𝑦 + 𝜐𝑥𝑦 + 𝑘𝑥𝑥 ) (3-37)
−ℎ/2 1 − 𝜐𝑥𝑦 𝜐𝑦𝑥
2 ∙ 𝐺𝑥𝑦 ∙ ℎ3
{ ∙ 𝑘𝑥𝑦
12
avec :
1 𝜕𝑢𝛼 𝜕𝑢𝛽 1 𝜕 2𝑤
𝜀𝛼𝛽 = ( + ) et 𝑘𝛼𝛽 = ( );
2 𝜕𝑥𝛽 𝜕𝑥𝛼 2 𝜕𝑥𝛼 𝜕𝑥𝛽
𝐸𝑥 et 𝐸𝑦 sont les modules d’Young (Pa) selon les directions x et y respectivement, 𝜐𝑥𝑦 et
𝜐𝑦𝑥 les coefficients de Poisson entre les directions 𝑥 et 𝑦, 𝐺𝑥𝑦 le coefficient de cisaillement
et ℎ l’épaisseur de la plaque (m).
𝜕 4 𝑤(𝑥, 𝑦) ′
𝜕 4 𝑤(𝑥, 𝑦) 𝜕 4 𝑤(𝑥, 𝑦)
𝐷𝑥 + 2(2𝐷𝑥𝑦 + 𝐷 𝑥𝑦 ) + 𝐷𝑦 + 𝑝(𝑥𝑚 ) = 0 (3-38)
𝜕𝑥 4 𝜕𝑥 2 𝜕𝑦 2 𝜕𝑦 4
avec :
𝐸𝑥 ∙ ℎ3 𝐸𝑦 ℎ3 𝐺𝑥𝑦 ∙ ℎ3 𝜐𝑦𝑥 ∙ 𝐸𝑥 ∙ ℎ3
𝐷𝑥 = ; 𝐷𝑦 = ; 𝐷𝑥𝑦 = ; 𝐷′𝑥𝑦 =
12(1 − 𝜐𝑥𝑦 𝜐𝑦𝑥 ) 12(1 − 𝜐𝑥𝑦 𝜐𝑦𝑥 ) 12 12(1 − 𝜐𝑥𝑦 𝜐𝑦𝑥 )
87
Chapitre III : Modélisation mécanique
E∙ℎ3
et 𝐷= la rigidité de flexion.
12(1−υ2 )
La méthode de Ritz peut être utilisée pour déterminer une solution approchée pour
des conditions aux limites plus générales, tant que des fonctions d’approximation
appropriées peuvent être trouvées pour le problème. La principale limitation de la
méthode de Ritz est que la solution converge lentement, à moins que la complexité des
fonctions d’approximation ne soit augmentée.
Une méthode analytique, Analytical Strip Method (ASM) [4], a été retenue. Cette
méthode analytique permet l’analyse de plaques soumises à différents types de
chargement (ponctuel, linéique, partiellement uniforme, surfacique, non uniforme). Elle
repose sur une décomposition de la plaque en bandes, dont le nombre dépend des
discontinuités géométriques, du type et de la localisation des chargements.
88
Chapitre III : Modélisation mécanique
La méthode repose sur une séparation des variables, où la solution pour chaque
bande de la plaque est exprimée par une décomposition en série de Fourier et les
charges sont exprimées par une série correspondante. Le principal avantage de la
méthode est que, contrairement aux méthodes numériques, la précision du résultat ne
dépend pas de la discrétisation du système, mais du nombre de modes considérés dans
la série de Fourier. Le nombre d’équations et la difficulté de calcul sont fortement réduits
grâce à cette solution semi-analytique.
Figure 3-12. Définitions des bandes et poutres dans le modèle mécanique plaque.
avec la solution homogène 𝑊ℎ𝑖 (𝑥, 𝑦) et 𝑊𝑝𝑖 (𝑥, 𝑦) la solution particulière qui dépend
directement du type de chargement appliqué à la bande 𝑖. Toutes les solutions possibles
pour 𝑊ℎ𝑖 (𝑥, 𝑦) sont données en ANNEXE E. On y trouve aussi les solutions particulières
𝑊𝑝𝑖 (𝑥, 𝑦) pour les chargements les plus courants. Si la bande 𝑖 est soumise à une
89
Chapitre III : Modélisation mécanique
𝑑 4 𝑊𝑝𝑗
𝑄𝑝𝑗 = 𝐸𝑝𝑗 𝐼𝑝𝑗 (3-42)
𝑑𝑦 4
𝑑2 θ𝑝𝑗 𝑑 4 θ𝑝𝑗
𝑇𝑏𝑗 = 𝐺𝑝𝑗 𝐽𝑝𝑗 ∙ + 𝐸𝑝𝑗 C𝑝𝑗 ∙ (3-43)
𝑑𝑦 2 𝑑𝑦 4
90
Chapitre III : Modélisation mécanique
𝜕𝑊𝑖 𝜕𝑊𝑖+1
= 𝜕𝑊𝑖 𝜕𝑊𝑖+1
𝜕𝑥 𝜕𝑥 θ𝑝𝑗 = − =−
𝜕𝑥 𝜕𝑥
𝑀𝑥𝑖 = −𝑀𝑥(𝑖+1) 𝑄𝑝𝑗 = 𝑅x(i+1) − 𝑅𝑥𝑖 + 𝐽𝑚𝑖
𝑅𝑥𝑖 = 𝑅𝑥(𝑖+1) 𝑇𝑝𝑗 = 𝑀𝑥(𝑖+1) − 𝑀𝑥𝑖
𝒙=𝟎 𝑊1 = 0 𝑊1 = 0 𝑀𝑥1 = 0
𝑊𝑝1 = 𝑊1
𝑀𝑥1 = 0 𝜕𝑊1 𝑅𝑥1 = − 𝐽𝑚0 𝜕𝑊1
=0 θ𝑝1 = −
𝜕𝑥 𝜕𝑥
𝑄𝑝1 = 𝑅𝑥1 + 𝐽𝑚0
𝑇𝑝1 = − 𝑀𝑥1
𝑊𝑁 = 0 𝑊𝑁 = 0 𝑀𝑥𝑛 = 0
𝑊𝑝𝑆 = 𝑊𝑁
𝒙 = 𝒙𝑵 𝜕𝑊𝑁
𝑀𝑥𝑁 = 0 𝑅𝑥𝑁 = −𝐽𝑚𝑛 𝜕𝑊𝑁
=0 θ𝑝𝑆 = −
𝜕𝑥 𝜕𝑥
𝑄𝑝𝑆 = −𝑅𝑥𝑁 + 𝐽𝑚𝑛
𝑇𝑝𝑆 = −𝑀𝑥𝑁
La superposition est employée quand une crête est soumise à une combinaison
du chargement afin de déterminer la fonction de chargement sur l’arête :
91
Chapitre III : Modélisation mécanique
Le modèle et son implémentation ont été soumis à différents cas test pour
validation pour ensuite adapter le modèle au cas du dispositif étudié.
𝒂 = 0,084 m
𝒃 = 0,223 m
𝒉 = 500 µm
x1 = a/5
x2 = 2a/5
x3 = 3a/5
x4 = 4a/5
𝑬𝒙 = 𝑬𝒚 = 71 GPa
𝝏𝒙 = 𝝏𝒚 = 0,3
𝑮 = 25,8 GPa
𝑳𝒚𝟏 = 1000/b N/m
𝑳𝒚𝟐 = 500/b N/m
𝑳𝒚𝟑 = 1500/b N/m
𝑳𝒚𝟒 = 200/b N/m
- calcul avec une discrétisation en bandes et quatre charges définies au niveau des
quatre arêtes (cas 1),
- application du principe de superposition sur les déformations obtenues par le
calcul des quatre sous-problèmes :
o chargement sur 𝑳𝒚𝟏 , 2 bandes et 4 appuis, (cas 2),
o chargement sur 𝑳𝒚𝟐 , 2 bandes et 4 appuis, (cas 3),
o chargement sur 𝑳𝒚𝟑 , 2 bandes et 4 appuis, (cas 4),
o chargement sur 𝑳𝒚𝟒 , 2 bandes et 4 appuis, (cas 5).
92
Chapitre III : Modélisation mécanique
On constate que les déformations obtenues avec ces deux méthodes, Figure 3-14
𝑏
pour une ligne en 𝑦 = 2, sont identiques. Par la suite, la méthode de calcul de la
déformation de la plaque dans son ensemble sera retenue, car moins couteuse en temps
de calcul.
Figure 3-14. Déformation au niveau de la ligne y = b/2 des 5 situations – application du théorème de
superposition.
93
Chapitre III : Modélisation mécanique
Déformation de la plaque
0.14
0.2
0.12
X: 0.042
Y: 0.1115
Z: 0.1588
0.15
0.1
W(m)
0.1
0.08
0.05
0.06
0
0 0.04
0.02
0.25
0.04 0.2 0.02
0.06 0.15
0.1
X(m) 0.08
0.05 Y(m)
0.1 0
Plaque Poutre
𝑎 = 0,084 m 𝑎 = 3 mm
𝑏 = 0,223 m 𝑏 = 0,223 m
ℎ = 500 µm ℎ = 75 µm
𝐸𝑥 =𝐸𝑦=71GPa 𝐸𝑐𝑢𝑖𝑣𝑟𝑒 =115 GPa
𝜕𝑥 = 𝜕𝑦 = 0,3
𝐺= 25,8 GPa
𝐿𝑦𝑛 = 5761 N/m
Figure 3-16. Données d’entrée pour l’évaluation du modèle plaque avec 4 chargements linéiques et 4
poutres en cuivre.
La plaque est soumise à quatre chargements linéiques, et des conditions aux limites de
type encastrement ont été définies sur les quatre bords de la plaque. Il est visible (Figure
3-17) que la prise en compte des conducteurs sur le calcul de la déformation de la plaque
ne modifie que très légèrement le déplacement maximal obtenu.
94
Chapitre III : Modélisation mécanique
0.04 0.04
0.035 0.035
0.03 0.03
X: 0.042 X: 0.042
0.05 Y: 0.1115
0.05 Y: 0.1115
0.025
Z: 0.04445 0.025 Z: 0.04422
W(m)
0.02 0.01 0.02 0.01
0 0
0 0.25
0.25 0
0.02 0.02 0.2
0.2
0.04 0.15 0.04 0.15
0.06 0.1 0.06 0.1
Y(m)
X(m) 0.08 0.05 Y(m) X(m) 0.08 0.05
0.1 0 0.1 0
Figure 3-17. Déformations issues du modèle plaque avec 4 chargements sans a) et avec b) des
poutres en cuivre et 4 encastrements comme conditions aux limites.
Figure 3-18. Position des plots de rappel, vue supérieure sur le plan XY.
Pour ce faire, chaque plot est assimilé à un ressort dont la raideur équivalente est
obtenue à l’aide des caractéristiques des plots. La force élastique résultante est ajoutée
dans la fonction de chargement 𝐽𝑚𝑖 (Equation 3-44) comme force ponctuelle.
95
Chapitre III : Modélisation mécanique
𝑆∙𝐸
avec 𝑘𝑒𝑞 =
𝑙
-3
x 10
0.22
4
0.2
0.18 3.5
0.16
3
0.14
2.5
Y(m)
0.12
0.1 2
0.08 1.5
0.06
1
0.04
0.5
0.02
Figure 3-19. Comparaison pour la prise en compte des plots par des ressorts équivalents : l’approche
analytique à gauche, modèle COMSOL® à droite.
Les résultats sont très proches, avec un déplacement maximal de 4,5 mm pour l’approche
éléments finis et de 4,4 mm avec l’approche analytique.
Les Figure 3-20 à Figure 3-22 donnent des détails de ces résultats le long des lignes :
96
Chapitre III : Modélisation mécanique
-3
x 10
4.5
3.5
Méthode analytique
Méthode FEM
3
2.5
W (m)
2
1.5
0.5
0
0 0.05 0.1 0.15 0.2 0.25
Y (m)
Figure 3-20. Comparaison approche éléments finis et semi-analytique sur une ligne au niveau des
conducteurs.
-3
x 10
3.5
2.5
Méthode analytique
Méthode FEM
2
W (m)
1.5
0.5
0
0 0.05 0.1 0.15 0.2 0.25
Y(m)
Figure 3-21. Comparaison approche éléments finis et semi-analytique sur une ligne au niveau des
plots.
97
Chapitre III : Modélisation mécanique
-3
x 10
4.5
3.5
3
Méthode analytique
Méthode FEM
2.5
W (m)
2
1.5
0.5
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
X (m)
Figure 3-22. Comparaison approche éléments finis et semi-analytique sur une ligne transversale au
milieu du dispositif.
Les résultats mettent en avant l’importance des plots dans la déformation de la plaque,
via l’action de rappel qu’ils exercent sur la plaque.
98
Chapitre III : Modélisation mécanique
-3
x 10
3.5
m=5
3 m=6
m=7
m=8
2.5 m=9
m=10
m=11
2 m=12
W(m)
m=13
m=14
1.5
0.5
0
0 0.05 0.1 0.15 0.2 0.25
Y(m)
Figure 3-23. Déplacement sur une ligne au niveau des plots en fonction du nombre de termes dans la
somme.
-3
x 10
5
m=5
4.5 m=6
m=7
4 m=8
m=9
m=10
3.5
m=11
m=12
3
m=13
W(m)
m=14
2.5
1.5
0.5
0
0 0.05 0.1 0.15 0.2 0.25
Y(m)
Figure 3-24. Déplacement sur une ligne au niveau des conducteurs en fonction du nombre de termes
dans la somme
450
400
Nombre d'inconnues
350
300
250
200
150
5 6 7 8 9 10 11 12 13 14
m
III.5 Conclusions
Dans un premier temps, étant donné qu’un motif élémentaire constitué d’un seul
conducteur (avec son vis-à-vis comme source d’effort) avait pu être mis en évidence lors
de la modélisation électromagnétique, c’est ce dernier qui a guidé la modélisation
mécanique vers une approche poutre. Les hypothèses simplificatrices de cette
modélisation, n’ont pas permis de converger vers des résultats exploitables en terme de
déplacement. Néanmoins, cette approche a permis de mettre en avant l’impact du
feuilletage du conducteur en deux lamelles.
Une approche de type plaque a ensuite été mise en œuvre. Une modélisation
analytique reposant sur un découpage en bandes a été retenue et progressivement
enrichie. Comparativement à une résolution éléments finis, il a été mis en avant que
l’approche analytique reproduit bien le comportement mécanique du tapis de dégivrage.
Seulement, l’amplitude des déplacements obtenue reste éloignée des mesures effectuées
par le partenaire Zodiac (ceci est également vrai avec la solution éléments finis). Il
convient cependant de noter que tous les calculs ont été effectués avec la valeur
maximale du courant de décharge (courant pic de 5800 A) en statique. Comme il sera
montré par la suite, la dynamique du système met en avant un fort contraste entre les
constantes de temps électrique et mécanique. Le différentiel entre les constantes de
100
Chapitre III : Modélisation mécanique
101
Chapitre III : Modélisation mécanique
102
Chapitre 4 :
Etude énergétique et modèle
électrique
Chapitre 4 : Etude énergétique et modèle électrique
IV.1 Introduction..............................................................................................................105
IV.3.e Analyse des pertes par courants de Foucault parasites au niveau de la surface
à dégivrer ....................................................................................................................111
IV.4 Assemblage des modèles et approche statique avec boucle de calcul ....................113
IV.1 Introduction
Le cahier des charges fait mention d’un niveau d’énergie de 0,9 𝐽/𝑐𝑚² nécessaire
pour obtenir un déplacement de 350 µ𝑚 [3]. Cette donnée a été calculée en prenant en
compte l’énergie totale fournie pendant les essais en entrée et ramenée à la surface
totale du tapis. En d’autres mots, c’est le rendement global du dispositif.
Les dimensions du tapis étudié sont de 89 𝑚𝑚 𝑥 231 𝑚𝑚, soit une surface de
205.6 𝑐𝑚². L’énergie théorique totale 𝑊 stockée dans les condensateurs d’alimentation
peut être calculée à partir de l’équation 4-1 :
1
𝑊 = 𝐶𝑉 2 (4-1)
2
105
Chapitre IV : Etude énergétique et modèle électrique
L’enjeu de cette partie est d’arriver à estimer la part des différentes énergies dans
le but de savoir si :
Le circuit électrique génère une impulsion de courant durant 1ms avec des crêtes
de 4 kA à 10 kA (de durée 40 µs à 50 µs), issue de la décharge de condensateurs de
capacité 800 µF-2000 µF préalablement chargé sous une tension de 400 VDC à 800 VDC
[5].
106
Chapitre IV : Etude énergétique et modèle électrique
Ainsi, il est possible de proposer le schéma électrique de la Figure 4-2 sur lequel
apparaissent les éléments parasites du serpentin, le modèle du condensateur ainsi que la
charge du système. Cette charge, modélisée à l’aide d’une impédance équivalente, traduit
l’énergie mécanique consommée nécessaire à la déformation.
La nature de cette impédance peut être exprimée par l’association d’un élément résistif en
série avec une inductance : la résistance représente les pertes, l’inductance représente la
dynamique du dispositif.
Dans une logique de conception d’un modèle viable pour les démarches
d’optimisation, une expression analytique des grandeurs électriques est recherchée.
107
Chapitre IV : Etude énergétique et modèle électrique
Le schéma électrique équivalent du système établi (Figure 4-2) conduit à une équation
différentielle du second ordre pour le courant :
𝑑2 𝑖(𝑡) 𝑑𝑖(𝑡)
(𝑙𝑠 + 𝐿𝑐ℎ𝑎𝑟𝑔𝑒)𝐶 + (𝑅𝑐ℎ𝑎𝑟𝑔𝑒 + 𝑟𝑠 + 𝑟𝑐)𝐶 + 𝑖(𝑡) = 0 (4-2)
𝑑𝑡 2 𝑑𝑡
Par intégration de l’équation 4-2, l’expression de la forme d’onde du courant peut être
obtenue. La résolution analytique est atteignable pour l’équation qui régit ce circuit, la
solution analytique de (Equation 4-2) pour ce circuit est définie par (Equation 4-3) :
𝑉0 − 𝑉𝑇
𝑖(𝑡) = (𝑒 𝛼1𝑡 − 𝑒 𝛼2𝑡 ) (4-3)
(𝑙𝑠 + 𝐿𝑐ℎ𝑎𝑟𝑔𝑒)(𝛼1 − 𝛼2)
1
𝛼1 = −
√(𝑙𝑠 + 𝐿𝑐ℎ𝑎𝑟𝑔𝑒)𝐶
2
−(𝑅𝑐ℎ𝑎𝑟𝑔𝑒 + 𝑟𝑠 + 𝑟𝑐)𝐶 − √(𝑅𝑐ℎ𝑎𝑟𝑔𝑒 + 𝑟𝑠 + 𝑟𝑐)𝐶 − 4(𝑙𝑠 + 𝐿𝑐ℎ𝑎𝑟𝑔𝑒)𝐶
𝛼2 =
2(𝑙𝑠 + 𝐿𝑐ℎ𝑎𝑟𝑔𝑒)𝐶
108
Chapitre IV : Etude énergétique et modèle électrique
Pour modéliser la charge, les valeurs des mesures disponibles (tension et courant
d’alimentation), issues des essais effectués par notre partenaire industriel, sont utilisées.
Pour compléter l’analyse énergétique, on s’appuie aussi des résultats de modélisation
précédents.
L’énergie absorbée par le tapis : 𝐸𝑎𝑏𝑠 = ∑ 𝑉(𝑖∆𝑡) ∗ 𝐼(𝑖∆𝑡) ∗ ∆𝑡 𝐸𝑎𝑏𝑠 = 283 J soit 1,38 J/cm²
𝑖
109
Chapitre IV : Etude énergétique et modèle électrique
Du Tableau 4-1, il ressort que pour une énergie absorbée de 283 J soit 1,38 J/cm²,
une énergie mécanique de 20,5 J est obtenue. Ramenée à la surface totale du tapis, cela
correspond à une densité de 0,1 J/cm². Donc, 0,1 J/cm² d’énergie mécanique ont été
suffisants pour assurer la fonction désirée du système, c’est-à-dire, la déformation de la
surface du tapis.
110
Chapitre IV : Etude énergétique et modèle électrique
Il est montré que la modélisation effectuée est satisfaisante. Les équations pour le modèle
électrique sont validées.
Pour compléter l’étude énergétique, il s’avère nécessaire d’évaluer les pertes par
courants de Foucault induites sur le blindage en aluminium. Ce dernier, étant à proximité
du dispositif de dégivrage, est affecté par le rayonnement du champ électromagnétique.
111
Chapitre IV : Etude énergétique et modèle électrique
Figure 4-5. Distribution des courants de Foucault [A/m2] induites sur la plaque pour le premier
harmonique à 625 Hz.
Le Tableau 4-2 liste les valeurs des pertes issues des courants induits dans la plaque
pour chaque harmonique.
Tableau 4-2. Pertes totales dans la plaque dues aux courants induits liées à l’impulsion de courant.
112
Chapitre IV : Etude énergétique et modèle électrique
A partir de cette décomposition harmonique des pertes par courant de Foucault, il est
possible d’estimer l’énergie associée sur un cycle de fonctionnement (Tableau 4-3).
On observe qu’au niveau de la zone à dégivrer, les pertes associées aux courants
de Foucault développées dans la plaque sont négligeables. Elles ne représentent donc
pas une grandeur pertinente à considérer en vue d’une optimisation de la consommation
du dispositif.
113
Chapitre IV : Etude énergétique et modèle électrique
Une première approche repose sur une succession d’états statiques dans une
boucle temporelle, selon le schéma proposé Figure 4-6. Boucle évaluée n fois, avec n le
nombre d’échantillons de la courbe de courant.
𝑛 =0; 𝑤 =0
𝑡𝑛+1 = 𝑡𝑛 +△ 𝑡
Calcul de l’amplitude du
courant
𝐼(𝑡𝑛+1 )
Evaluation de l’effort
électromagnétique (calcul CADES)
𝐹(𝑡𝑛+1 )
𝑤(𝑡𝑛+1 )
Elle a été mise en œuvre sur la géométrie initiale du dispositif. Les résultats
obtenus sont présentés dans la Figure 4-7 et montrent que cette approche permet de
retrouver les premiers résultats pour la valeur du déplacement.
114
Chapitre IV : Etude énergétique et modèle électrique
Avec cette approche, il est observé que la dynamique de la déformation suit celle
du courant. Ceci est conforme à l’utilisation d’une approche statique, ne prenant pas en
compte la dynamique de la réponse mécanique pour le déplacement. Ainsi, on observe
que cette approche ne reflète pas totalement la réalité puisque le retard de la réponse
mécanique (l’inertie mécanique), par rapport à la sollicitation électrique, n’apparait pas. Il
s’en suit une sur évaluation du déplacement obtenu.
-3
x 10
1.2
Point A
Point B
1 Point C
Point D
Déplacements(m)
0.8 Point E
0.6
0.4
0.2
0
0 1 2 3 4 5 6
temps(s) x 10
-4
Figure 4-7. Evolution du déplacement en fonction du temps pour les 5 points de référence sur le tapis
(d’après la Figure 1-8).
Le modèle électrique conçu peut être couplé au reste des modèles précédemment
développés.
115
Chapitre IV : Etude énergétique et modèle électrique
A ce stade, l’origine de ces pertes ne peut pas être identifiée avec plus de
précision mais elle pourrait se trouver au niveau des pertes Joule dans les conducteurs,
ou encore, être liées aux restrictions mécaniques issues de la constitution structurelle du
tapis dégivreur.
Les valeurs des éléments du circuit électrique équivalent ont été obtenues à l’aide
des mesures. La principale limitation reste sur l’identification des éléments équivalents de
la plaque qui ne tient pas compte du développement des conducteurs lors d’un cycle
d’activation.
116
Chapitre IV : Etude énergétique et modèle électrique
117
Chapitre 5 :
Analyse dynamique et résultats
Chapitre 5 : Analyse dynamique et résultats
V.2 Synthèse des analyses de sensibilité réalisées au cours du projet sur la base des
modèles numériques développés....................................................................................122
V.3.b Essai expérimental 2 : le déplacement est imposé par les forces génères au
niveau des lignes des conducteurs ..............................................................................139
V.5 Conclusions..........................................................................................................145
Chapitre V : Analyse dynamique et résultats
V.1 Introduction
Dans un premier temps, des études de sensibilité, menées en utilisant les modèles
précédemment établis, sont présentés. Par la suite une analyse des données
expérimentales fournies par Zodiac va permettre de visualiser l’évolution des grandeurs
en fonction du temps et cerner les détails du comportement du déplacement en fonction
du temps.
121
Chapitre V : Analyse dynamique et résultats
Des études de sensibilité ont été lancées au cours du projet avec les différents
modèles indépendamment : électromagnétique, mécanique, électrique. Le but étant
d’avoir des premières analyses de sensibilité.
Il s’avère utile d’identifier les enjeux qui peuvent avoir de l’importance dans la
définition des grandeurs optimales de conception pour le dispositif. Une typologie des
paramètres de l’espace de recherche d’optimisation est proposée :
Paramètres topologiques
- Le nombre de lignes de conducteurs,
- la densité de plots,
- définition du type de jonction aux extrémités du tapis,
- disposition du circuit serpentin.
Paramètres géométriques
- Les dimensions géométriques du tapis,
- les dimensions géométriques des conducteurs et des plots,
- l’épaisseur du blindage en aluminium.
Paramètres physiques
- Caractéristiques mécaniques, tels que le module de Young ou le module de
cisaillement, des matériaux des conducteurs, du blindage, des élastomères.
Paramètres d’excitation
- Forme d’onde du courant : amplitude, temps de montée, temps de descente,
- stratégie de commande : une seule impulsion, plusieurs impulsions, etc.
122
Chapitre V : Analyse dynamique et résultats
Une question légitime est de savoir comment faire évoluer la géométrie du tapis
pour à la fois optimiser la consommation énergétique et s’assurer de fournir l’énergie
demandé, de 0.1 J/cm², pour la transformation mécanique.
D’un point de vue topologique, pour une surface totale du tapis fixe, il est possible de faire
varier :
Le but est de trouver les paramètres influents en vue d’identifier des pistes pour
l’optimisation du dispositif.
Dans tous les cas, on cherche comment va évoluer la déformation du tapis, ainsi
que les circuits équivalents et la force électromagnétique résultante.
Pour le modèle électrique utilisé dans les tests suivants, une forme d’onde de
courant avec un modèle, ne prenant pas en compte l’élément impédance de charge
couplé, est utilisé.
Cette étude préliminaire a été effectuée avant la définition d’un modèle électrique
final mais elle est toutefois utile dans les objectifs d’analyse du fonctionnement. Elle
permet en effet d’apprécier les variations directement liées aux modifications
géométriques imaginées pour le dispositif.
123
Chapitre V : Analyse dynamique et résultats
Deux conducteurs supplémentaires sont insérés, comme illustré sur la Figure 5-1.
Par ailleurs, pour une même surface de tapis, dans cette configuration, le nombre de
lignes de chargement mécanique est augmenté, ainsi que le nombre de lignes de plots de
rappel.
a) b)
Une modélisation InCa3D a permis d’évaluer les nouvelles valeurs du schéma électrique
équivalent du dispositif : 𝑅𝑠𝑒𝑟𝑝𝑒𝑛𝑡𝑖𝑛 = 105 𝑚𝛺 et 𝐿𝑠𝑒𝑟𝑝𝑒𝑛𝑡𝑖𝑛 = 138 𝑛𝐻, valeurs considérées
constantes pour la gamme de fréquences étudiée.
124
Chapitre V : Analyse dynamique et résultats
a) b)
125
Chapitre V : Analyse dynamique et résultats
Le Tableau 5-2 présente les principaux résultats obtenus pour cette configuration.
Tableau 5-2. Comparaison entre le cas initial et la structure avec conducteurs de largeur double.
La Figure 5-4 présente le déplacement des différents points pris le long d’une ligne
au milieu du dispositif (Figure 5-1 b) en fonction du temps. Du fait du courant plus élevé,
la force électromagnétique est plus grande et donc le déplacement plus important. Cette
configuration est donc très intéressante. L’inconvénient majeur réside dans
l’augmentation de la masse de cuivre nécessaire, qui serait intéressant d’évaluer.
126
Chapitre V : Analyse dynamique et résultats
Cette configuration est intéressante mais il faudra veiller à ne pas trop diminuer le
nombre de plots dont la présence a un intérêt pour la structure, qui est de garantir le
retour des conducteurs à leur position initiale après chaque activation du dégivreur.
Tableau 5-3. Comparaison entre le cas initial et la structure avec moins de plots.
127
Chapitre V : Analyse dynamique et résultats
Tableau 5-4. Comparaison entre le cas initial et la structure avec un tapis moins large.
Dans ce cas, le déplacement maximal (Figure 5-6) a diminué car le fait de réduire la
surface du tapis rend la structure plus rigide. Certes, par rapport au tapis original les
pertes ont diminué mais l’augmentation de la densité d’énergie surfacique n’a pas été
suffisante pour contrecarrer l’augmentation de la rigidité. Aussi, la valeur de la densité
d’énergie surfacique n’est pas le seul indicateur de la possible adéquation du tapis réalisé
au cahier des charges.
Figure 5-6. Déformation dans le cas d’un tapis moins large (cas 4).
128
Chapitre V : Analyse dynamique et résultats
maximal. Le test mise en œuvre propose de vérifier si la souplesse gagnée compense les
pertes énergétiques pour une longueur de conducteurs augmentée de 50 %. Le Tableau
5-5 montre les principaux résultats obtenus.
Tableau 5-5. Comparaison entre le cas initial et la structure avec des conducteurs plus longs.
Tout d’abord le nombre de plots est resté constant (cas 5, Figure 5-7). Dans ce
cas, le tapis est plus souple, non seulement parce qu’il est rallongé mais aussi car il y a
moins de plot rapportés à la surface. Le cas 5bis représente la même configuration mais
avec une densité de plots constante (Figure 5-8).
Figure 5-7. Déformation dans le cas de conducteurs plus longs (cas 5).
129
Chapitre V : Analyse dynamique et résultats
Figure 5-8. Déformation dans le cas de conducteurs plus longs (cas 5bis).
130
Chapitre V : Analyse dynamique et résultats
Il est observé une influence importante sur les valeurs de déplacement pour le
changement de la taille du tapis et pour le nombre de plots. La géométrie des
conducteurs reste aussi un paramètre de poids parce qu’elle est directement liée aux
grandeurs électriques telles que la résistance et en conséquence le courant.
Pour cette étude, le modèle mécanique présenté dans le paragraphe III.4.d, est
pris en compte. Les calculs sont effectués pour une densité de force linéique fixe (𝐿𝑦𝑖 =
5400 N/m) et des conditions aux limites de type encastrement sur les 4 bandes de la
plaque sont considérées.
Pour une géométrie et un matériau donnés, le nombre de plots par ligne varie de 3
à 10, les plots sont disposés de manière uniforme sur 3 lignes parallèles pour un nombre
de termes m dans la somme de l’équation 3-34, suffisamment élevé pour s’affranchir
d’une quelconque erreur liée à ce facteur a été pris en compte (m=14). Sur la Figure 5-9,
les résultats pour le déplacement le long d’une ligne au niveau d’un conducteur sont
représentés.
131
Chapitre V : Analyse dynamique et résultats
-3
x 10
16
3 plots
14 4 plots
5 plots
6 plots
12
7 plots
8 plots
10 9 plots
10 plots
8
W(m)
-2
0 0.05 0.1 0.15 0.2 0.25
Y(m)
Figure 5-9. Déformation le long d’une ligne au niveau d’un conducteur en fonction du nombre de plots.
Il est constaté que plus le nombre de plots est important, plus le déplacement est faible.
La force de rappel exercée par les plots limite le déplacement. Le déplacement maximal
obtenu en fonction du nombre de plots confirme ce résultat (Figure 5-10)
0.016
0.014
0.012
Wmax (m)
0.01
0.008
0.006
0.004
3 4 5 6 7 8 9 10
Nombre de plots
132
Chapitre V : Analyse dynamique et résultats
Cependant, si le nombre de plots est trop réduit, la force de rappel permettant au tapis de
retrouver sa forme d’origine ne sera pas suffisante. Cela pourrait impacter la fiabilité du
dispositif après plusieurs cycles d’activation.
Figure 5-11. Relation entre le nombre de plots, la densité de force et le courant d’alimentation
garantissant une valeur de déplacement imposée ( 𝑾𝒎𝒂𝒙 = 𝟎, 𝟎𝟎𝟒𝟒 𝒎 ).
Une relation de proportionnalité est clairement apparente : pour assurer une valeur de
déplacement Wmax donnée, le courant d’alimentation doit être d’autant plus élevé que le
nombre de plots est grand. La relation directe entre le courant d’alimentation et la densité
de force s’explique facilement du fait de l’équation de force de Laplace (Equation 2-3).
133
Chapitre V : Analyse dynamique et résultats
Une comparaison a été effectuée une nouvelle fois avec COMSOL® pour
confirmer la fiabilité du comportement du modèle développé. La Figure 5-12 montre
clairement l’impact des plots sur la déformation du tapis, avec la présence d’une surface
ondulée.
Figure 5-12. Comparaison pour la prise en compte des plots comme des ressorts : approche
analytique, modèle Comsol à droite.
134
Chapitre V : Analyse dynamique et résultats
-3
x 10
4.5
3.5
3
Wmax (m)
2.5
1.5
0.5
0 0.5 1 1.5 2 2.5 3 3.5 4
Module d'Young du matériau élastomère (Pa) x 10
7
Figure 5-13. Déplacement maximal en fonction du module d’Young du matériau du plot élastomère.
0.015
0.01
Wmax (m)
0.005
0
1 2 3 4 5 6 7 8 9
Diamètre des plots (m) x 10
-3
135
Chapitre V : Analyse dynamique et résultats
Ces études de sensibilité ont permis d’éliminer certaines variables, soit parce
qu’elles se sont avérées avoir peu d’influence soit parce qu’elles sont déjà fixées à leur
valeur optimale ou à cause des contraintes de fabrication. Seule une optimisation à
travers les deux modèles résultants (modèle électrique et modèle électro-magnétique-
mécanique), permettra de faire une définition plus précise des grandeurs optimales pour
les variables déjà identifiées.
Il a été observé que l’amplitude du courant est un paramètre important : plus celle-
ci est élevée, plus l’effort sera important et donc le déplacement également, moyennant
de bien prendre en compte les pertes.
De la même façon, l’écart entre conducteurs est un paramètre de grand impact sur
le déplacement.
Il est à noter que les différentes analyses de sensibilité sur les paramètres du tapis
ont été effectuées avec un modèle d’impédance de charge 𝑍𝑐ℎ𝑎𝑟𝑔𝑒 constant, déduit de
l’identification à partir des mesures (et sans prise en compte de la charge 𝑍𝑐ℎ𝑎𝑟𝑔𝑒
équivalente dans le cas des études de sensibilité du paragraphe V.2.a). Les paramétres
des conducteurs sont pris en compte dans le modèle PEEC (𝐿𝑠𝑒𝑟𝑝𝑒𝑛𝑡𝑖𝑛 , 𝑅𝑠𝑒𝑟𝑝𝑒𝑛𝑡𝑖𝑛 ) pour
toutes les études de sensibilité présentées. Cependant, à ce stade l’impact de la
déformation du tapis sur l’impédance de charge, qui traduirait le couplage
électromagnétique, n’est pas pris en compte.
136
Chapitre V : Analyse dynamique et résultats
A partir des résultats des expériences réalisées avant la thèse, et dans le but de
mieux cerner les caractéristiques du comportement dynamique de la réponse mécanique
et sa relation avec le signal électrique d’alimentation, deux séries d’essais effectuées par
le partenaire Zodiac sont reprises ici.
Ces deux essais ont été faits à vide, c’est-à-dire en l’absence de glace, afin de pouvoir
mesurer le déplacement maximal fourni par le dispositif. Les résultats complets issus de
ces essais sont disponibles dans les documents de référence et rapports, propriétés de
Zodiac Aerospace.
Compte tenu des imprécisions dans la procédure de mesure, pour cette première
série d’essais, seul le comportement de la déformation est analysé. Les valeurs de
déplacement obtenu ne sont pas considérées pour le moment.
Sur la Figure 5-15 et Figure 5-16, le déplacement d’une ligne au droit des plots et
des conducteurs a été tracé en fonction du temps.
137
Chapitre V : Analyse dynamique et résultats
Figure 5-15. Déplacement d’une ligne au droit des plots en fonction du temps.
Figure 5-16. Déplacement d’une ligne au droit d’un conducteur en fonction du temps.
- à l’exception des extrémités les déplacements sont uniformes, tous les points
de la ligne se déplacent de la même quantité. Temporellement cela se produit
pendant la phase de croissance du courant. Il s’agit donc d’un régime forcé de
déformation,
- par la suite, les points de la ligne se déplacent de façon plus aléatoire. Le courant
décroit jusqu’à s’annuler. Il s’agit d’un régime libre au cours duquel la structure
libère l’énergie emmagasinée.
138
Chapitre V : Analyse dynamique et résultats
Une autre série d’essais a été effectuée de manière plus précautionneuse sur 5
points d’une même ligne transversale (Figure 5-17 et Figure 5-18). Les points A, C et E se
situent au niveau de lignes de plots alors que les points B et D sont au niveau des
conducteurs internes du dispositif. Pour cet essai, une série de 5 décharges identiques de
deux condensateurs en parallèle chargés sous 600 VDC a été réalisée.
A
B
C
D
E
Figure 5-18. Position des points de mesure sur le dispositif de dégivrage (conducteurs et plots), plan
XY du système.
Les Figure 5-19 et Figure 5-20 présentent les courbes de déplacement obtenues
pour les 5 points en fonction du temps. Au niveau des conducteurs, une certaine symétrie
est observée. Au niveau des plots, cette symétrie est moins marquée pour deux des trois
points, A et E (points aux extrémités). Le point C qui est au milieu du tapis a un
déplacement maximal moindre comparé à celui des points B et D au niveau des lignes de
conducteurs, ce qui est cohérent compte tenue de la conception du tapis.
139
Chapitre V : Analyse dynamique et résultats
Il est constaté une plus grande cohérence et plus d’uniformité dans les courbes
relevées à partir de ces essais. Pour cette raison, les résultats de cet essai serviront de
référence afin de valider les hypothèses émises et alimenter les analyses.
140
Chapitre V : Analyse dynamique et résultats
liée au courant de décharge qui augmente progressivement n'est pas suffisante pour
contrer la résistance mécanique imposée par la plaque, les conducteurs ne se déplacent
pas.
𝜕 2 𝑤(𝑥, 𝑦, 𝑡)
𝐷 ∙ Δ2 𝑤(𝑥, 𝑦, 𝑡) + 𝛾ℎ − 𝑝(𝑥, 𝑦, 𝑡) = 0 (5-1)
𝜕𝑡 2
La solution analytique n’étant pas abordable, la stratégie retenue est la résolution par
différences finies.
Dans la stratégie abordée jusqu’à présent, la boucle de calcul statique, très facile
à mettre en œuvre, ne permet pas de résoudre l’équation 5.1. Pour pallier ce défaut, une
autre stratégie est d’utiliser une résolution numérique basée sur la méthode des
différences finies avec un maillage spatial bidimensionnel et un algorithme de Newmark
dû à la présence d’une équation différentielle d’ordre 2 en temps.
142
Chapitre V : Analyse dynamique et résultats
Figure 5-21. Algorithme de résolution par la méthode aux différences finies [2].
Cette procédure a été mise en œuvre sur le dispositif Zodiac étudié et les résultats
obtenus pour les déplacements des mêmes 5 points (A, B, C, D, et E) sont présentés sur
la Figure 5-22.
143
Chapitre V : Analyse dynamique et résultats
Figure 5-22. Evolution du déplacement obtenu par la méthode aux différences finies.
Les courbes montrent l’apparition d’un retard au démarrage lié au terme temporel de
l’équation 5.1. Comparativement aux résultats expérimentaux, l’allure des courbes
obtenues au début du cycle est satisfaisante. Il existe une symétrie de la réponse entre
les points A-E et B-C, ces courbes se superposent.
144
Chapitre V : Analyse dynamique et résultats
Figure 5-23. Evolution du déplacement avec la méthode par différences finies et la force de contact.
Une augmentation de la valeur du déplacement maximal est constatée. Ceci peut être dû
au changement dans la discrétisation.
V.5 Conclusions
145
Chapitre V : Analyse dynamique et résultats
146
Chapitre V : Analyse dynamique et résultats
147
Conclusion générale et
perspectives
Conclusions et perspectives
152
Conclusions et perspectives
151
Conclusions et perspectives
des conducteurs. Les résultats qualitatifs ont permis de conclure que les hypothèses
simplificatrices applicables et validées pour la modélisation électromagnétique ne sont pas
toutes transférables pour la description mécanique. Une approche de modélisation
analytique par déformation de plaques avec une discrétisation en bandes, a ensuite été
développée. Des résultats plus satisfaisants ont conduit à sélectionner cette approche et
continuer les travaux de définition du modèle; celui-ci a été enrichi avec la définition
d’éléments structurels du tapis. Ces calculs restent cependant dans un cadre statique.
Comme premier pas vers l’optimisation, dans le chapitre 5 une série d’études de
sensibilité est présentée. Parmi les paramètres topologiques sensibles identifiés, on trouve
la section des conducteurs, le nombre de lignes de charge, la densité et diamètre des plots
de rappel. Des études de sensibilité concernant les paramètres dits d’extraction ont aussi
152
Conclusions et perspectives
montré le grand impact de l’écart entre les conducteurs et du courant d’alimentation sur les
résultats.
Le besoin d’une prise en compte de la dynamique est mise en évidence avec une
étude sur le comportement des formes d’onde de déformation mesurées lors des essais
réalisés précédemment à la thèse. Pour répondre à ce besoin, une modélisation complète
basée sur une résolution par la méthode des différences finies a été implémentée. Elle
permet la prise en compte de la variation temporelle de la sollicitation électrique et de la
réponse mécanique. Grâce à cette procédure, une réponse dynamique plus proche des
essais a été obtenue. Pour compléter ces travaux de modélisation, une contrainte de
contact pour limiter le déplacement en cohérence avec la structure du dispositif a été
ajoutée.
L’étude réalisée dans ces travaux de thèse permet d’affirmer la pertinence d’une
technologie électromécanique pour le dégivrage aéronautique. Malgré une modélisation
multiphysique complexe, cette solution permet d’estimer les caractéristiques de
déformation obtenue pour l’action de dégivrage. Si une optimisation est mise en œuvre,
l’efficacité énergétique de cette solution électromécanique pourrait aussi s’avérer meilleure
que celle des solutions traditionnelles, du fait que l’énergie est directement transformée et
appliquée à la surface à traiter. Néanmoins, il reste à évaluer les possibilités en termes de
solutions d’alimentation adaptées à la fois aux contraintes aéronautiques et aux fortes
intensités de courant demandées par cette solution.
153
Conclusions et perspectives
Perspectives
154
Conclusions et perspectives
155
Annexes
Annexes
i(
On continu par étudier l’influence du nombre des fils du serpentin sur l’induction
électromagnétique B, paramètre directement proportionnel à la déformation de notre
surface à dégivrer.
La Figure A-2 présente les premiers trois cas à analyser, avec la définition des
paramètres d (distance entre conducteurs extérieures selon l’axe Y) et e (écart sur l’axe Z
entre les deux serpentins). Dans un premier instant, on augmente le nombre de fils en
respectant toujours la taille totale du tapis (donc les fils s’approchent quand leur nombre
augmente), dans la deuxième analyse, on augmente le nombre de fils mais on va faire
varier la taille du tapis (la distance entre conducteurs reste égale).
159
Annexes
Coupe
X=>cte
Figure A-2. Augmentation du nombre de fils conducteurs pour une même valeur d.
Avant de passer aux équations, il est nécessaire d’établir la convention des axes,
directions et variables à utiliser, ainsi que la géométrie du système et le sens des
courants.
Tableau A-1. Description des données et variables d’étude.
Si l’on se repère sur la Figure 2-5, on souhaite analyser la force qui va générer un
déplacement dans la direction Z. De ce fait c’est la composante Y du champ magnétique B
d’intérêt. Remarquer que la contribution By de n’importe quel élément situé à la même
position Z sera nulle (θ=90°, cos θ=0) ; l’induction électromagnétique By sur un fil
conducteur sera donc celle créée par la contribution des champs de tous les fils placés
dans une position Z différente au fil analysé, en autres termes la force subie par le
serpentin du haut n’est due qu’au champ B créé par les conducteurs du serpentin du bas
160
Annexes
Situation #1
d y
µ I µ I 3 µ I
cos y
0 0 0
By 1( y ) arctan cos arctan cos
2 2 e 2 e 2
2 y e d y e2 2 d y e2
2 2
3 3
d y 2 d y
3 µ I 3 µ I
d y
cos
0 0
s arctan cos arctan arctan
e 2 e 2 2 e
2 2 d y e2 2 ( d y ) e
3
Situation #2
d y
µ I µ I 4 µ I
cos y
0 0 0
By 2( y ) arctan cos arctan cos ar
2 2 e 2 e 2
2 y e d y e2 2 d y e2
2 2
4 4
d y 2 d y 3 d y
4 µ I
0 4 µ I
0 4 µ I
0
s arctan cos arctan cos arctan
e 2 e 2 e 2
2 2 d y e2 3 d y e2 2 ( d y )
2
4 4
3 d y
4 µ I
ctan
0
cos d y
arctan
e 2 2 e
2 ( d y ) e
Situation #3
d y
µ I µ I 5 µ I
cos y
0 0 0
By 3( y ) arctan cos arctan cos arct
2 2 e 2 e 2
2 y e d y e2 2 d y e2
2 2
5 5
d y 2 d y 3 d y
5 µ I
0 5 µ I
0 5 µ I
0
os arctan cos arctan cos arctan c
e 2 e 2 e 2
2 2 d y e2 3 d y e2 4 d y e2
2 2
5 5 5
161
Annexes
3 d y 4 d y
5 µ I 5 µ I
d y
cos arctan
0 0
n cos arctan
e 2 e 2 2 e
2 4 d y e2 2 ( d y ) e
5
Figure A-3. Courbes résultantes pour les trois situations étudiées (By1 Induction magnétique pour le
cas de 4 fils, By2 pour le cas de 5 fils, By3 pour le cas au 6 fils).
Il peut être remarqué que le fait d’augmenter le nombre de fils n’assure pas une
modification de l’induction magnétique générée. En effet, l’induction magnétique est trop
faible et s’annule entre les conducteurs, ce qui tend à penser que l’influence du
conducteur le plus proche est prépondérante.
162
Annexes
Sur la Figure A-4, on observe l’induction au niveau d’un seul conducteur isolé, celui au
bord du tapis dans les 3 cas précédents et on rajoute le premier terme de ces expressions.
On confirme bien que c’est uniquement le conducteur le plus proche qui contribue de
manière significative à l’induction. On peut alors négliger non seulement les conducteurs
qui se trouvent à la même position sur l’axe Z mais aussi tous les conducteurs qui
appartiennent à un autre tour du serpentin.
Figure A-4. Comparaison induction magnétique pour tout le système avec la contribution de la spire
plus proche.
163
Annexes
l’impédance des différentes configurations d’une part afin que l’on puisse savoir si la
demande énergétique varie avec les paramètres géométriques,
l’induction magnétique pour les situations 1, 2 et 3 pour faire la comparaison avec les
calculs analytiques.
On va utiliser InCa3D sur les géométries précédentes afin d’évaluer l’écart entre le cas
théorique (conducteurs filiformes de longueur infinie) et la géométrie réelle.
La Figure A-6 présente une superposition des cinq courbes résultantes pour la valeur de
l’impédance des différentes situations en prenant la configuration initiale comme référence.
164
Annexes
Figure A-6. Valeurs de l’impédance pour les cinq géométries analysées en fonction de la fréquence.
Prenons comme exemple la Figure A-7 qui montre une coupe au milieu du serpentin,
les iso valeurs montrées appartiennent au calcul de la densité de force de Laplace. On peut
clairement remarquer que cette distribution n’est pas uniforme sur toute la section, en
conséquence l’induction magnétique ne l’est pas non plus.
165
Annexes
Effectivement, dans le cas d’une section non nulle avec une densité de courant uniforme il
est aisé d’exprimer l’induction pour un conducteur de longueur infinie.
166
Annexes
Figure A-8. Edition du scénario de résolution, possibilité de piloter les paramètres géométriques établis
comme variables pendant la création de la géométrie.
167
Annexes
Figure A-10. Exemple de l’évaluation de la Densité de Force volumique sur une grille avec une courbe 3D.
Composant 3 => Axe Z de notre système.
Figure A-11. Possibilité de passer aux courbes 2D pour plus de simplicité et de clarté dans la
présentation des résultats.
168
Annexes
Grâce aux résultats des études précédentes, où on démontrait l’impact prépondérant des
conducteurs plus proches pour les calculs électromagnétiques, on a pris en compte
uniquement une section du système originalement construit par ZODIAC : deux conducteurs
longitudinaux alimentés par le même courant en sens opposé. En vue de simplifier la
modélisation mécanique vers l’analyse d’une poutre 2D, on a commencé à considérer un
seul conducteur d’épaisseur 150µm par couche au lieu de deux conducteurs d’épaisseur
75µm par couche soudés aux extrémités du dispositif original.
Cas d’analyse 1vs5 : Deux conducteurs parallèles, le premier constitué d’un seul élément et le
deuxième divisé en 5 éléments
169
Annexes
Conducteur_0
Figure A-12. Géométrie simulé pour le cas #2 : a) InCa3D vue 3D et plan XZ b) McMmems vue 3D.
Sur le Tableau A-2 on trouve les données d’entrée du premier scenario pour la géométrie du
cas 1vs5. Les résultats sont montrés sur les graphiques de Figure A-13 et Figure A-14. La
force par conducteur en fonction de la variation du courant DC d’alimentation du circuit est
présentée, ainsi que le pourcentage d’erreur estimé entre les deux approches de simulation.
170
Annexes
Tableau A-2. Valeurs d’entrée pour les paramètres géométriques et physiques des simulations (Cas
1vs5, 1er scenario).
171
Annexes
Figure A-13. Force de Laplace produite sur les conducteurs en fonction du courant d'alimentation (Cas #2. 1vs5, 1er scenario).
Annexes
Figure A-14. Pourcentage d’erreur entre les calculs des deux logiciels (Cas #2. 1vs5, 1er scenario).
Les risques sont la manipulation de plusieurs fichiers pour passer d’un outil de la
plateforme à l’autre. La création des dossiers n’est pas automatisée et peut favoriser les
erreurs humaines.
Pour plus de clarté, la nouvelle séquence de simulation peut être décrite suivant les
étapes énoncées ci-dessous (Tableau A-3):
- Création d’un fichier XML (extension *.xml) avec le logiciel McMmems. Dans ce
document on détermine la géométrie et les paramètres du système à évaluer ainsi
que des équations additionnelles si on le souhaite.
- Depuis la même interface McMmems, on exporte le modèle vers un fichier SML
(extension *.sml), code contenant une série d’équations auto générée.
173
Annexes
Logo et
ordre de 1 2 Densité
de Force
chainage
Plate-
McMMems CADES Generator CADES Calculator Matlab
forme
Exten-
sion du *.sml *.icar
*.xml *.m
compo- *.icar *.val
sant
Interface pour
l’édition du code qui Interface pour la
Logiciel pour le
Interface pour la définit la géométrie, définition des valeurs
codage des
définition de la ses paramètres, d’entrée des
fonctions de calcul
géométrie, ses unités, équations et paramètres
mathématique
caractéristiques variables de sortie. précédemment
diverses. A l’aide
Descrip- physiques et Sa compilation définis, lancer des
d’un plug-in dédié,
tion énonciation des origine les calculs, et tracer des
il est possible de
calculs requis. Une composants ICAR courbes. Les
piloter les
fois compilé, on qui pourront ensuite données d’entrée
composants ICAR
génère le composant être lus par la sont sauvegardées
et récupérer ses
SML. calculette CADES en forme de fichier
sorties.
ou par un VAL.
programme Matlab.
174
Annexes
175
Annexes
176
Annexes
177
Annexes
178
Annexes
Les coefficients Cm1, Cm2, Cm3, Cm4, détaillés dans les Tables A-4 et A-5, dépendent du
type de conditions aux limites sur les deux arêtes (y = 0 et y = b) et de la valeur de m
(nombre de termes de la somme infinie). Ces coefficients, la fonction Gm(y) et g2m sont
déterminés à partir du Tableau A-6.
179
Annexes
180
Annexes
Les coefficients Cm1, Cm2, Cm3, Cm4 et la fonction Gm(y) sont aussi déterminés à partir du
Tableau A-6.
181
Annexes
Le cas d’une plaque chargée ponctuellement a été traité ici. Dans un premier temps, ce ne
sont pas les valeurs exactes des déplacements qui seront observées mais plutôt le
comportement de ce modèle et la vérification de certains cas particuliers pour valider la
cohérence des résultats.
a = 0,084 m
b = 0,223 m
Epaisseur = 5x10-4 m
Ex = Ey = 71x109 Pa
µx= µy = 0,3
G = 25,8x109 Pa
P = 1000 N
182
Annexes
Figure A-15. Déplacement de la plaque soumise à une charge ponctuelle au niveau des lignes (x= a/2 et
y = b/2), 4 appuis.
Déformation de la plaque
0.12
X: 0.042
Y: 0.1115
Z: 0.14
0.14
0.1
0.12
0.1
0.08
0.08
W(m)
0.06 0.06
0.04
0.02 0.04
0
0
0.02 0.25 0.02
0.2
0.04 0.15
0.06 0.1
X(m) 0.08 0.05 Y(m)
0.1 0
Figure A-16. Déformation de la plaque en 3D soumise à une charge ponctuelle de 1000 N, 4 appuis.
183
Annexes
On dispose de deux références bibliographiques [1] et [2], pour lesquelles ce cas a été
traité, le premier à l’aide d’une approche éléments finis utilisant l’outil numérique Abaqus
et le deuxième article avec une approche analytique. On s’est servi de ces deux
références pour faire une comparaison avec notre modèle et vérifier que l’implémentation
a été bien réalisée. Le Tableau A-9 résume les résultats pour ce cas test traité avec les
trois approches différentes.
Tableau A-9. Résultats pour le déplacement maximal de la plaque chargée ponctuellement sur le centre
résolu avec 3 méthodes différentes.
Notre approche analytique (avec m =5) FEM [6] Autre approche analytique [7]
Afin de vérifier la symétrie de l’approche, on a refait les calculs avec les mêmes données
d’entre mais maintenant pour une plaque carrée (a = b = 0,223 m). On prend encore une
fois comme référence le tracé pour la déformation de la plaque sur deux lignes centrales.
Les résultats de ce cas sont présentés sur la Figure A-17.
0.8
Déformation de la ligne y = b/2
en fonction de x
0.7 Déformation de la ligne x = a/2
en fonction de y
0.6
0.5
W(m)
0.4
0.3
0.2
0.1
Figure A-17. Déformation sur deux lignes centrales d’une plaque carrée (a = b = 0,223 m), charge
ponctuelle de 1000 N, 4 appuis.
Les deux courbes sont identiques ce qui met en évidence la symétrie de la déformation
entre les deux axes X et Y. On conclut alors que les différences rencontrées entre les
courbes de déformation trouveés pour la plaque rectangulaire de la Figure A-15 sont en
184
Annexes
effet liées à la différence entre les mesures a≠b, et non à une erreur de représentation du
modèle. De cette façon, on valide ce premier modèle et son implémentation.
La notion de poutre est nécessaire dans ce modèle car elle sert à définir l’endroit où un
effort linéique est appliqué. Pour notre dégivreur, les poutres seront indispensables pour la
définition des efforts linéiques génères au niveau des lignes de conducteurs du serpentin.
Ce cas, issu d’une référence bibliographique [3] a été sélectionné pour caler
l’implémentation des éléments poutre. Cette fois, sur deux côtés opposés de la plaque
carré, on définit deux poutres pour remplacer les conditions aux limites.
CC 0,4322 0,4344
CS 0,5400 0,5446
SS 0,7106 0,7098
Effectivement, notre modélisation converge avec les résultats de la littérature, ce qui veut
dire que les caractéristiques physiques des poutres ont bien été prises en compte.
On utilise le premier modèle de plaque rectangulaire avec les mêmes données d’entrée du
Tableau A-8. Sur ce cas simple, on veut maintenant regarder l’impact du découpage en
bandes au niveau de la nouvelle arête ainsi définie.
185
Annexes
convergence des résultat par rapport au premier l’analyse pour la plaque rectangulaire non
discrétisée de la Figure A-15. L’implémentation de la discrétisation est alors validée.
Figure A-18. Déformation au niveau des lignes (x= a/2 et y = b/2) pour une plaque avec une charge
ponctuelle au centre de 1000 N, discrétisation en 2 et 4 bandes, CL= 4 appuis.
Avec :
Ici la valeur est fixée à 𝑚 = 5. Chaque bande a 4 inconnues : AmI, BmI, CmI et DmI (voir
Table A3 de l’ANNEXE F). Pour les cas de notre exemple, le nombre d’inconnues est :
On rappelle que le dégrée de précision du résultat pour cette méthode ne dépend pas du
nombre d’éléments de discrétisation ou bandes mais du nombre de termes de la série
considères 𝑚 . On a alors intérêt à limiter le nombre de bandes de discrétisation au
minimum nécessaire pour la définition de toutes les charges et autres caractéristiques de
notre système.
On continue avec le modèle plaque rectangulaire avec charge centrale ponctuelle, sauf
que pour ce cas de test les conditions aux limites de la plaque seront modifiées. L’objectif
de cette étude est d’évaluer la relation entre le type d’ancrage et la déformation de la
plaque.
186
Annexes
Les deux cas extrêmes identifiés pour encadrer le résultat sont l’appui simple et
l’encastrement. Sur les Figure A-19 à Figure A-23, les résultats de la déformation sont
montrés pour différentes positions sur la plaque.
0.1
W(m)
0.08
0.06
0.04
0.02
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
X(m)
Figure A-19. Déformation de la plaque pour différentes conditions aux limites pour les deux axes
d’analyse y=b/2 et x=a/2.
187
Annexes
Figure A-20. Déformation de la plaque en 3D du cas : Charge ponctuelle P = 1000 N, 2 appuis (sur deux
arêtes y = 0 et y = b) + 2 encastrements (sur deux arêtes x = 0 et x = a)
188
Annexes
Figure A-22. Déformation sur des lignes proches des bords, charge ponctuelle P = 1000 N, 2 appuis
(arêtes y = 0 et y = b) + 2 encastrements (arêtes x = 0 et x = a).
-4 -7
x 10 x 10
3 3
ligne y = 0
ligne x = 0
ligne y = b/5000
ligne x = a/5000
2.5 2.5 ligne y = b/4000
ligne x = a/4000
ligne y = b/3000
ligne x = a/3000
ligne y = b/2000
ligne x = a/2000
2 2
1.5 1.5
W(m)
W(m)
1 1
0.5
0.5
0
0
-0.5
-0.5 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
0 0.05 0.1 0.15 0.2 0.25 X(m)
Y(m)
Figure A-23. Déformation des lignes proches des bords, charge ponctuelle P = 1000 N, 2
encastrements (sur deux arêtes y = 0 et y = b) + 2 appuis (sur deux arêtes x = 0 et x = a).
Les images courbes 3D sont très parlantes. Pour tous le cas, le changement des
conditions aux limites a bien été pris en compte. On peut conclure que pour le modèle
actuel on sait contrôler et définir les conditions aux limites.
Afin de vérifier qu’il n’y a aucune erreur masquée par la symétrie du cas test traité, cette
fois on va faire varier la position du chargement ponctuel sur la plaque pour tester si la
189
Annexes
forme de la réponse est cohérente. Les deux cas test du Tableau A-11 ont été évalués
pour le modèle de plaque rectangulaire.
Tableau A-11. Changement de la position du chargement ponctuel sur la plaque.
Cas 1 Cas 2
On observe sur les Figure A-24 et Figure A-25, la forme de la surface de réponse pour les
déformations engendrées dans les deux cas. Le sens physique est bien conservé, ce qui
démontre la bonne utilisation de cette approche.
Déformation de la plaque
0.11
0.1
0.12 0.09
0.08
0.1
0.07
0.08 0.06
W(m)
0.06 0.05
0.04
0.04
0.03
0.02 0.02
0 0.01
0
0.02
0.25
0.04 0.2
0.06 0.15
X(m) 0.08 0.1
0.05 Y(m)
0.1 0
Figure A-24. Déformation de la plaque en 3D pour une charge ponctuelle de 1000 N au point (x = a/2, y
= b/5), 4 appuis
190
Annexes
Déformation de la plaque
0.06
0.05
0.08
0.04
0.06
0.04 0.03
W(m)
0.02 0.02
0
0.01
-0.02
0.25 0
0.2 0.1
0.15 0.08
0.06
0.1
Y(m) 0.04
0.05 0.02
X(m)
0 0
Figure A-25. Déformation de la plaque en 3D pour une charge ponctuelle de 1000 N au point (x = a/5, y
= b/5), 4 appuis
Maintenant que le modèle de plaque rectangulaire crée a passé tous les tests de
vérification avec un chargement concentré sur un point, il est temps d’augmenter la
complexité de la définition de la charge, vers une charge de type linéique. Pour ce faire, on
va se servir de la définition d’éléments poutre car c’est à travers eux qu’on peut
représenter les efforts mécaniques appliqués sur une ligne.
Dans cette configuration, Figure A-26, le type de chargement est modifié. On conserve la
même force totale de 1000 N pour cette occasion elle est répartie de façon uniforme le
long d’une ligne centrée, Ly=1000 N/m.
a = 0,084 m, b = 0,223 m,
Epaisseur = 5E-4 m
Ex = Ey = 71E9 Pa
µx= µy = 0,3
G = 25,8E9 Pa
Ly = 1000/b N /m
191
Annexes
La déformation a été évaluée le long des deux lignes x = a/2 et y = b/2 tel que pour la
batterie de cas de test 1, ces résultats sont présentés sur les courbes de la Figure A-27.
Les tapis de déformation 3D des figures Figure A-28 et Figure A-29 correspondent aux
résultats pour les analyses avec des conditions aux limites de type appuis simples ou
encastrement respectivement.
Figure A-27. Déformation de la plaque soumise à une charge uniformément répartie sur la ligne x=a/2,
au niveau des lignes centrées
Déformation de la plaque
0.06
X: 0.042 0.05
0.07 Y: 0.1115
Z: 0.06193
0.06
0.04
0.05
0.04
W(m)
0.03
0.03
0.02
0.02
0.01
0 0.25
0 0.2 0.01
0.02 0.15
0.04
0.06 0.1
0.08 0.05
X(m) 0.1 0
Y(m)
Figure A-28. Déformation de la plaque pour une charge de 1000/b N/m sur la ligne x = a/2, 4 appuis.
192
Annexes
Déformation de la plaque x 10
-3
16
14
0.02
X: 0.042 12
Y: 0.1115
Z: 0.01739
0.015 10
8
W(m)
0.01
6
0.005
4
0 2
0
0.02 0.25
0.04 0.2
0.06 0.15
X(m) 0.1
0.08 0.05 Y(m)
0.1 0
Figure A-29. Déformation de la plaque pour une charge de 1000/b N /m sur la ligne x = a/2, 4
encastrements.
Afin d’analyser ceci, une charge de 1000N uniformément répartie est appliquée,
mais cette fois uniquement sur un segment de la longueur totale 2∆L sur la ligne x = a/2
comme le montre la Figure A-30. Quand ∆L tend vers zéro on tend vers la configuration
d’un chargement ponctuel et quand ∆L tend vers b/2 on se rapproche du chargement sur
une ligne complète.
Figure A-30. Chargement de la plaque sur une ligne partielle : définition des grandeurs
193
Annexes
S’il y a convergence de cette approche avec ces deux cas extrêmes de charge ponctuelle
et charge linéique antérieurement traités on peut alors imaginer n’utiliser que le cas de
chargement partiel sur une ligne pour définir tous les cas.
0.08
W(m)
0.06
0.04
0.02
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
X(m)
Figure A-31. Déformation au niveau de la ligne y = b/2 : Charge Ly = 1000/(2∆L) N/m sur le morceau 2∆L
de la ligne x = a/2
Les résultats du test montrés sur la Figure A-31 confirment la continuité obtenue grâce à
cette formulation unique d’une charge uniforme sur une ligne partielle.
Dans ce paragraphe, tout comme dans le cas du chargement ponctuel l’ajout de poutres
pour remplacer deux conditions aux limites a été effectué. Dans [5], le cas présenté a été
traité et les valeurs maximales des déplacements sont donc comparées sur le tableau ci-
après. Notre modélisation converge avec les résultats de [5].
CL Wmax(m) Ref.
CC 0,0328 0,0334
CS 0,0426 0,0423
SS 0,0572 0,0559
Figure A-32. Plaque carrée, chargement sur une ligne partielle, 2 poutres.
194
Annexes
Dans cette partie, la position du chargement est modifiée pour supprimer la symétrie du
problème et éliminer d’éventuelles erreurs d’application de la méthode. Les charges sont
appliquées uniformément sur la ligne complète aux trois positions différentes : x = a/4, x =
a/2 et x = 3a/4. Les résultats sont présentés sur les Figure A-33 et Figure A-34 pour une
évaluation de la déformation sur des lignes centrées pour différents types de chargement.
Les Figure A-35 et Figure A-36, présentent la déformation en 3D.
Figure A-33. Déformation au niveau de la ligne y = b/2 avec 3 positions de charge sur une ligne
complète : x = a/4, x= a/2 et 3a/4.
195
Annexes
Figure A-34. Déformation au niveau de la ligne x = a/2 avec 3 positions de charge sur une ligne
complète : x = a/4, x= a/2 et 3a/4.
Déformation de la plaque
0.045
0.04
0.04
0.035
X: 0.042
0.035
0.03 Y: 0.1115
Z: 0.0424
0.025
W(m)
0.03
0.02
0.015
0.025
0.01
0.005
0.02
0
0.25
0.015
0.2
0.15 0.01
0.1
Y(m) 0.005
0.05 0.09
0.07 0.08
0.05 0.06
0.03 0.04
0 0.02
0 0.01
X(m)
196
Annexes
Déformation de la plaque x 10
-3
0.01
8
0.008
7
0.006
W(m)
X: 0.042
Y: 0.1115
0.004 Z: 0.008673 6
0.002
5
0
0.25 4
0.2 3
0.15
2
Y(m) 0.1
1
0.1
0.05 0.08
0.06
0.04
0 0.02
0
X(m)
Pour ce cas d’évaluation, on va changer les conditions aux limites qui sont définies sur les
arêtes au bord de la plaque, tel que sur le cas 1, à la différence que cette fois on ne vas
pas appliquer une charge ponctuelle mais la charge linéique. Les résultats de leur
influence sont montrés sur les Figure A-37 à la Figure 40.
Figure A-37. Déformation au niveau de la ligne y = b/2 en fonction des conditions aux limites de la
plaque, Ly = 1000/b (N/m) sur x = a/2.
197
Annexes
Figure A-38. Déformation au niveau de la ligne x = a/2 en changeant les conditions aux limites de la
plaque, Ly = 1000/b (N/m) sur x = a/2.
Déformation de la plaque
0.055
X: 0.042
0.06 Y: 0.1115
0.05
Z: 0.05737
0.05 0.045
0.04
0.04
0.035
W(m)
0.03
0.03
0.02 0.025
0.02
0.01
0.015
0
0 0.01
0.25
0.05 0.2 0.005
0.15
X(m) 0.1
0.05
0.1 0 Y(m)
Figure A-39. Déformation de la plaque : 2 encastrements (y=0 et y=b) et 2 appuis (x=0 et x=a), Ly =
1000/b (N/m) sur x = a/2.
198
Annexes
Déformation de la plaque x 10
-3
16
0.02
14
X: 0.042
0.015 Y: 0.1115
Z: 0.01733
12
10
W(m)
0.01
8
0.005
6
0 4
0.4
0.3 0.1
0.08 2
0.2 0.06
0.1 0.04
Y(m) 0.02
0 0
X(m)
Figure 40. Déformation de la plaque : 2 appuis (y=0 et y=b) et 2 encastrements (x=0 et x=a), Ly = 1000/b
(N/m) sur x = a/2.
Les résultats confirment que, pour ce cas aussi, le changement des conditions aux limites
influence fortement la déformation de la plaque. Le type de sujétion définie pour les
extrémités de notre système reste un paramètre fortement sensible. Paramètre que
d’ailleurs on ne connaît pas avec certitude.
Dans cette partie, l’influence d’un chargement non uniforme sera évaluée. Pour y parvenir
il faut une définition par segments sur la longueur de la ligne.
199
Annexes
Figure A-41. Définition du par segments pour un chargement linéique non uniforme.
Si les valeurs des chargements uniformes s’égalisent, les résultats convergent vers le cas
d’une seule charge sur la ligne complète. On peut ainsi faire une comparaison des
résultats avec ceux du cas de charge uniforme sur toute la ligne et valider cette nouvelle
implementation.
Il faut noter que quand l’arête xi est soumise à une combinaison de charges, la
méthode de superposition est utilisée pour déterminer la fonction 𝐽𝑚𝑖 .
𝒑
𝑱𝒎𝒊 (𝒚) = ∑𝒄=𝟏 𝑱𝒎𝒊𝒄 (𝒚) (A-2)
Où :
Les trois tests ci-dessous sont évalués avec les valeurs d’entré suivants :
Les résultats obtenus sont présentés sur les Figure A-42 à Figure A-45 pour les 3 cas en
2D et en 3D. On vérifie que le cas 1 tend vers le cas 3 quand les valeurs des charges non
uniformes s’homogénéisent ; comme espéré.
200
Annexes
Figure A-42. Déformation au niveau de la ligne y = b/2 pour les 3 cas considérés.
Figure A-43. Déformation au niveau de la ligne x = a/2 pour les 3 cas considérés.
201
Annexes
Déformation de la plaque
0.09
0.08
0.1
X: 0.042
0.07
Y: 0.1115
W(m)
Z: 0.08733
0.05
0.06
0.05
0
0 0.04
0.02
0.03
0.04
0.02
X(m) 0.06 0.25
0.2 0.01
0.08 0.15
0.1
0.05 Y(m)
0.1 0
Déformation de la plaque
0.06
0.07 X: 0.042
Y: 0.1115
Z: 0.06193 0.05
0.06
0.05
0.04
0.04
W(m)
0.03
0.03
0.02
0.01
0.02
0
0
0.25
0.02 0.2 0.01
0.04 0.15
0.06 0.1
X(m) 0.08 0.05 Y(m)
0.1 0
A ce stade on sait définir tous les types d’éléments dont on a besoin pour commencer à
définir la mécanique du système qu’on souhaite évaluer. On peut maintenant procéder à
l’adaptation de ce modèle pour l’application technologique du dégivreur électromécanique.
202
Annexes
Calcul de la FFT :
Ts 0.0016
1
Fs
Ts
Ne 300
Ts
T e
Ne
Ne
k 0 1
2
1
Fe
Te
Ne 1
I( n Te) exp i 2 k n
D( k ) Te
Ne
n 0
harmonique 1 harmonique 33
D( harmonique ) 0.969 3
D( harmonique ) 2.066 10
203
Annexes
Pour résoudre l’équation 3-39, une stratégie est d’utiliser une résolution numérique basée
sur la méthode des différences finies avec un maillage spatial bidimensionnel et un
algorithme de Newmark dû à la présence d’une équation différentielle d’ordre 2 en temps
(ou un algorithme Euler implicite pour une équation différentielle du premier ordre).
La première étape est de transformer l’équitation. 3-39 et pour cela les variables x,
y et t sont redéfinies comme suit, faisant ainsi apparaître les compteurs de boucles
spatiales (i et j) et temporelle (k) : x = ix ; y = jy ; t = kt.
𝜕4𝑤
( )
𝜕𝑥 2 𝜕𝑦 2 𝑖,𝑗,𝑘
En remplaçant toutes ces définitions dans Eq. 3-39 on obtient Eq. A-1 :
𝐷∆𝑡 2
𝑤𝑖,𝑗,𝑘+1 = − 𝜌
(𝑃𝑤𝑖,𝑗,𝑘 + 𝑄(𝑤𝑖+1,𝑗,𝑘 + 𝑤𝑖−1,𝑗,𝑘 ) + 𝑅(𝑤𝑖,𝑗+1,𝑘 + 𝑤𝑖,𝑗−1,𝑘 ) +
avec :
6 8 6 4 4 4 4
𝑃= 4
+ 2 2+ 4 ; 𝑄=− 4− 2 2 ; 𝑅=− 4− 2 2 ;
∆𝑥 ∆𝑥 ∆𝑦 ∆𝑦 ∆𝑥 ∆𝑥 ∆𝑦 ∆𝑦 ∆𝑥 ∆𝑦
2 1 1
𝑆 = ∆𝑥 2 ∆𝑦2 ; 𝑇 = ∆𝑥 4 ; 𝑈 = ∆𝑦4
𝑤|𝑥=0,𝑎 = 0 ; 𝑤|𝑦=0,𝑏 = 0
𝜕𝑤 𝜕𝑤
| =0; | =0
𝜕𝑥 𝑥=0,𝑎 𝜕𝑦 𝑦=0,𝑏
204
Annexes
Dans un premier temps, ce processus a été testé sur un exemple trouvé dans la littérature
afin de valider son implantation [4].
L’exemple traité est présenté Figure A-46 [8]. Il s’agit d’une plaque rectangulaire
avec un chargement ponctuel q.
Figure A-47, les résultats obtenus sont comparés à ceux de l’exemple trouvé dans la
référence bibliographique mentionnée. Il y a des différences en raison du choix de la
discrétisation dans le temps et dans l’espace, ce qui semble avoir une influence sur la
valeur maximale. Pour gérer la discrétisation en temps et dans l’espace on joue sur les
paramètres ∆𝑥, ∆𝑦 et ∆𝑡. Il faut aussi respecter la condition de convergence suivante ainsi
que le nombre de pas de temps qui gère la durée de la simulation (kmax=10000) :
1 𝐷∆𝑡 2 ∆𝑦 2 ∆𝑥 2
0≤𝑐≤ 𝑎𝑣𝑒𝑐 𝑐 = ( ) (2 + + )
4 𝜌∆𝑥 2 ∆𝑦 2 ∆𝑥 2 ∆𝑦 2
205
Annexes
Figure A-47. Comparaison sur un exemple : à gauche la littérature, à droite notre calcul.
Par la suite dans cet article, une FFT est réalisée. La comparaison des résultats est
présentée Figure A-48.
Figure A-48. Comparaison sur un exemple : FFT : à gauche la littérature, à droite notre calcul.
Des différences subsistent dues au fait que la durée de la simulation n’est pas suffisante.
Mais cet exemple permet tout de même de valider la démarche.
Pour finir il a été mis en œuvre sur le dispositif Zodiac étudié et les résultats
obtenus sont présentés Figure A-49. Le retard au démarrage lié au terme temporel dans
l’équation 3-39 est bien pris en compte et l’allure des courbes obtenues au début du cycle
sont satisfaisantes. Par la suite, lorsque le temps augmente des valeurs négatives de
déplacement apparaissent, ce qui montre que certaines contraintes n’ont pas été prises en
compte. Son analyse se poursuit.
206
Annexes
207
Annexes
Comme le montre par ailleurs la Figure A-51 ci-après, les conducteurs doivent
nécessairement se retrouver sur une ligne de discrétisation du maillage du modèle
différences finies : changer la position des conducteurs va nécessiter de modifier le
maillage.
208
Annexes
Figure A-51. Index pour construction du maillage en relation à la position des conducteurs et des plots.
Cela amené à développer une formulation avec un maillage déformable, dite approche par
différences finies à pas de discrétisation variable. Le modèle représente un préalable
nécessaire pour aller vers une optimisation se montrant quantitativement et
qualitativement plus réaliste. Un programme Matlab a été réalisé dans ce but.
Un des derniers verrous à lever en vue d’aller vers une optimisation du dispositif,
est de pouvoir s’appuyer sur un schéma électrique équivalent, dans lequel le tapis est
modélisé par une charge équivalente dont les paramètres dépendent des paramètres de la
structure en général, et du tapis en particulier.
209
Annexes
[1] Amar Khennane, Introduction to Finite Elements Analysis using Matlab and Abaqus.
CRC Press, 2013.
[2] S. Timoshenko and S. Woinowsky-Krieger, Theory of plates and shells. McGraw Hill,
1959.
[3] Issam E. Harik and Ghassan L. Salamoun, “The analytical strip method of solution for
stiffened rectangular plates,” Computers & Strucutres, vol. 29, no. 2, pp. 288–291,
1988.
[4] A.R. Tavakolpour, I. Z. M. Darus, and M. Mailah, “Numerical Simulation of a Flexible
Plate System for Vibration Control,” WSEAS TRANSACTIONS on SYSTEMS and
CONTROL, vol. 4, no. 5, Mar-2009.
210
FR - Modélisation et analyse d’un système de dégivrage en aéronautique
Ce manuscrit de thèse aborde les travaux de modélisation d’une nouvelle technologie électromécanique pour le
dégivrage en aéronautique en mode vol.
Le principe de fonctionnement de la technologie étudiée repose sur la déformation mécanique des surfaces à dégivrer
à partir d’efforts électromagnétiques générés par une excitation en courant électrique de forte intensité. La solution
imaginée par le partenaire industriel Zodiac Aerospace, a conduit à un démonstrateur dont les premiers essais ont
été lancés préalablement à cette thèse. À partir des résultats obtenus, l’objectif de ces travaux de thèse a été d’obtenir
un modèle adapté pour les futures démarches d’optimisation dans le but ultime d’aller vers une solution de dégivrage
plus électrique que les solutions existantes, performante, et plus efficace en termes de rendement énergétique.
Une modélisation multiphysique a été réalisée. Cette modélisation est constitué de plusieurs sous modèles
analytiques qui ont été intégrés sur une plateforme de résolution numérique. La modélisation mécanique met en
œuvre une approche de type plaque par éléments bande, qui a été résolu en fonction du temps par une méthode
de différences finies. Le modèle mécanique résultant est dit semi-analytique. La modélisation électromagnétique
repose sur une définition analytique et a présenté comme difficulté la complexité de la définition géométrique du
démonstrateur. Un circuit électrique équivalent du dispositif complet a été identifié et les équations du modèle
électrique ont été définies ; ceci a permis de réaliser une analyse énergétique pour comprendre la transformation des
forces de nature différente et mettre en évidence la possibilité d’une future optimisation.
Les résultats de simulation à partir du modèle final représentent assez bien la dynamique de la déformation
mécanique observée expérimentalement mais présentent en revanche des limitations en termes de précision.
This thesis report explains the modeling procedure of a new electromechanical de-icing technology in aeronautics for
in-flight application.
The operation principle of the technology resides in the mechanical deformation of the working surface caused by the
effect of electromagnetic forces generated from a high-intensity current source.
The industrial partner Zodiac Aerospace conceived this solution they carried out the construction of a demonstrator
and a series of tests prior to this Ph.D. Based on the obtained results, the aim of this project was to achieve an
adequate model of the de-icing system that should be suitable for further optimization of the device, such that a more
electric de-icing solution will be proposed, with good performance and with higher energetic efficiency.
A multiphysiscs model was developed, which comprises multiple analytical submodels that where integrated into a
numerical resolution platform. The mechanical submodel implements a strip approach for plates solved via a finite
differences method that permits time dependence evaluation, and can be defined as semi analytical. Another
submodel is based on the mathematical definition of the electromagnetic behavior, the main complication of which
was to consider the complex geometrical definition of the demonstrator. An equivalent electric circuit for the whole
system was identified and the equations for an electrical submodel where then established. This allowed the study of
the energy transformation and repartition, and reveals the possibility of future optimization.
Simulation results from the final model properly reproduce the dynamics of the mechanical deformation response as
were observed during the previous experiments, but also reveals some minor accuracy problems.