CHARLES VILLEMURE

OPTIMISATION À L'AIDE D'ALGORITHMES

GÉNÉTIQUES D'UN STRATIFIÉ POREUX SOUMIS
À UN FLUX THERMIQUE EN CONVECTION
NATURELLE

Mémoire présenté
à la Faculté des études supérieures de l'Université Laval
dans le cadre du programme de maîtrise en génie mécanique
pour l'obtention du grade de maître es sciences (M.Se.)

FACULTÉ DES SCIENCES ET DE GÉNIE

DÉPARTEMENT DE GÉNIE MÉCANIQUE
UNIVERSITÉ LAVAL
QUÉBEC

2007

© Charles Villemure, 2007

 Résumé

Cette étude traite de l'optimisation d'un échangeur de chaleur constitué d'une jux­
taposition de différentes couches de matériaux poreux. Grâce à un réseau interne de
pores, un fluide refroidissant est libre de circuler à travers l'échangeur par convection
naturelle. La problématique consiste à déterminer la distribution optimale de porosité
et l'ordonnancement idéal des matériaux de la série de couches qui minimisent la tempé­
rature maximale dans le refroidisseur. Au cours de cette analyse, la chaleur transmise au
système émane d'une plaque chaude adjacente au stratifié. Le flux thermique est donc
fixe sur cette frontière. Le champ de température et l'écoulement, déterminés selon
les caractéristiques physiques et dimensionnelles du stratifié, sont calculés de manière
numérique. Les équations de conservation de quantité de mouvement et d'énergie sont
résolues par une méthode basée sur les volumes finis. Le processus d'optimisation se
rattachant à la distribution de porosité et à l'assignation des matériaux dans le stratifié
est réalisé à l'aide d'un algorithme génétique (AG). De plus, une composition opti­
male du système, correspondant à un minimum global, peut être atteinte sous certaines
contraintes, s'exprimant notamment en matière de coût et de masse. Cette approche
permet également de dimensionner l'échangeur, car l'AG a la capacité d'éliminer une
portion du stratifié afin de répondre aux exigences demandées.

i
 Abstract

In order to meet modem engineering needs in terms of cooling Systems (e.g., cooling
of electronics), porous structures (e.g. metallic foams) are regarded as an interesting
alternative to fins due to their large surface of heat transfer per volume. In this paper,
we investigate the optimal configuration of a porous médium structure in order to reduce
its thermal résistance. The System studied represents a stacking of porous layers, in
which a cooling fluid circulâtes, adjacent to a heat-generating surface. The flow within
the staking is driven by natural convection. The objective of this work is to minimize
the hot spot température of the system. The design variables are the porosities and
materials of each layer. The thermal performance is evaluated with a CFD code based
on the finite volume approach. The hot spot température minimization is pursued under
global mass and cost constraints, with a genetic algorithm (GA). The GA détermines
the optimal porosities and solid material of each layer. Purthermore, the optimal total
length of the stacking is indirectly determined by the GA as layers can be added or
removed in order to improve the global performance and satisfy the constraints. The
results of the présent study reveal that an appropriate distribution of porosity and
material benefits the minimization of the température in a layered porous médium.

ii
"Un expert est une personne qui a fait toutes les erreurs qui peuvent être
faites dans un domaine très étroit"
Niels Bohr, 1885-1962

m
 Avant-propos

Tout d'abords, je tiens à remercier de façon toute particulière mon directeur de

recherche, Louis Gosselin, dont la disponibilité, la rigueur scientifique et la créativité ont
significativement contribué à mener à terme cette maîtrise. Je remercie également Guy
Gendron, mon co-directeur de recherche, pour les judicieux conseils qu'il m'a transmis
tout au long de ce travail.

Je souligne par ailleurs la contribution de mes collègues de travail. Je pense entre

autres à Guillaume, Mario, Marc, Yanick et Philippe avec qui la vie au laboratoire ne
fut que plus agréable.

Je tiens tout particulièrement à remercier ma famille pour leur support. Je suis

également reconnaissant envers Geneviève, dont l'amour et l'appui ont joué un rôle
incontestable dans l'accomplissement de ce projet.

IV
 Table des matières

Résumé i

Abstract ii

Avant-propos iv

Table des matières v

Liste des t a b l e a u x viii

Table des figures x

Liste des s y m b o l e s xiii

1 Introduction 1
1.1 Problématique 1
1.2 Revue bibliographique 3
1.3 Objectifs 6
1.3.1 Objectif principal 6
1.3.2 Objectifs secondaires 7
1.4 Méthodologie et présentation du document 7

2 Modèle mathématique 8
2.1 Présentation de la géométrie du système 9
2.2 Description du modèle mathématique 10
2.3 Solution analytique du régime couche limite 17
2.3.1 Approximation par ordre de grandeur du régime couche limite . 18

v
 2.3.2 Approche par la méthode intégrale du régime couche limite . . . 20
2.3.3 Porosité optimale en régime couche limite 24
2.4 Plage d'opération du modèle 24
2.5 Propriétés des matériaux employés 26
2.6 Conclusion 27

M o d è l e numérique 28
3.1 Maillage utilisé 29
3.2 Discrétisation par volumes finis 30
3.2.1 Exemple de résultats numériques 33
3.3 Indépendance de maillage 34
3.4 Validation du modèle numérique 37
3.5 Conclusion 44

A l g o r i t h m e génétique 45
4.1 Introduction 46
4.2 Mécanismes de l'algorithme génétique 48
4.2.1 Génération de la population initiale 48
4.2.2 Décodage 49
4.2.3 Correction des individus 51
4.2.4 Évaluation des individus 53
4.2.5 Classement des individus 54
4.2.6 Sélection 56
4.2.7 Reproduction 57
4.2.8 Recherche locale 58
4.2.9 Mutation 61
4.2.10 Nouvelle génération 63
4.2.11 Critère d'arrêt et convergence 64
4.3 Conclusion 65

R é s u l t a t s de l'optimisation 66
5.1 Optimisation sans contraintes 67
5.1.1 Influence de la longueur sur la résistance thermique 67
5.1.2 Effet du nombre de couches 72
5.1.3 Distribution de porosité 72
5.1.4 Masse et coût 75
5.2 Optimisation sous une contrainte de masse 77
5.3 Optimisation sous une contrainte de coût 81

vi
 5.4 Optimisation sous une combinaison de contraintes de masse et de coût 83

Conclusion 93

Bibliographie 96

A E x e m p l e de p r o b l è m e dimensionnel 101

B Script de l'algorithme génétique 103

B.l Résumé du contenu de l'AG 103
B.2 Boucle principale 105
B.3 Fonctions de l'AG 115
B.3.1 Banque de matériaux 115
B.3.2 Historique 116
B.3.3 Correspondance binaire-réel 118
B.3.4 Correspondance réel-binaire 120
B.3.5 Création de la population initiale 122
B.3.6 Correction des individus 124
B.3.7 Tronquer la séquence 126
B.3.8 Résolution thermofluide 128
B.3.9 Conductivité équivalente 138
B.3.10 Perméabilité 140
B.3.11 Résolution d'un système triagonal 141
B.3.12 Intégration d'une fonction 142
B.3.13 Pénalité de la fonction objectif 143
B.3.14 Classement 145
B.3.15 Recherche locale 147
B.3.16 Croisement multi-points 154
B.3.17 Mutation 156
B.3.18 Réinsertion 157
B.3.19 Sélection des individus 160

vii
 Liste des tableaux

2.1 Propriétés des matériaux disponibles lors de l'optimisation 26

3.1 Etude d'indépendance de maillage pour un stratifé fait d'aluminium,

L = l,(j) = 0.7, N = 1, D = 0.001 35
3.2 Étude d'indépendance de maillage pour un stratifé fait d'aluminium,
l = 0.5, 0 = 0.95, N =1, D = 0.001 36
3.3 Ny requis pour l'indépendance de maillage en fonction du Ra et de 0
pour un stratifé fait d'aluminium, L = 0.5, N = 1, D = 0.001 36
3.4 Longeur requis pour un développement couche limite en fonction du Ra,
0 = 0.5 et 0 = 0.9 , N = 1,D = 0.001 38
b
3.5 Coefficients de l'expression f (y) = ay , 0 = 0.5, iV = 1, L = 5,L> = 0.001 42
3.6 Coefficients de l'expression Tw(y) = ayb , <p = 0.90, N = 1, L = 5,
D = 0.001 42
3.7 Validation du modèle numérique en considération de T m a x ,0 = 0.5, N =-
1, D = 0.001 42
3.8 Validation du modèle numérique en considération de T m a x , 0 = 0.90,
N = 1, D = 0.001 43
3.9 Validation du modèle numérique en considération de Tmoy, 0 = 0.90,
N = 1, L = 5, D = 0.001 43
3.10 Validation du modèle numérique en considération de Tmoy, 0 = 0.5,
jV = 1, L = 5, £> = 0.001 44

5.1 Caractéristique et structure du design optimal pour : Ra — 10 16 ,L =

0.5,iV = 10, em = ec = 0, D = 0.001 71

vm
5.2 Longeur efficace en fonction du nombre de Rayleigh, N = 10, em = ec ==
0, D = 0.001 71

IX
 Table des figures

1 Schéma de l'échangeur 10
2 Représentation schématique d'un empilement de micro-conduits . . . . 12
3 Développement de la couche limite thermique et température maximale
dans le stratifié 20

1 Exemple de maillage : L = 1, N = 4, Ny = 16 29
2 Discrétisation par volumes finis 30
3 Distribution de température et profil de vitesse dans un échangeur po­
reux, tf> = 0.5, Ra = 10 15 , L = 5, D = 0.001 33
4 Distribution de température et profil vitesse dans un échangeur poreux,
(f) = 0.5, Ra = 10 15 , L = 5, D = 0.001 34
5 Exemple du maillage utilisé avec le logiciel Fluent : L = 1, N = 1,
Nv = 16, D = 0.001 39
15
6 Validation du modèle : T en fonction de y, <f> — 0.5, Ra = 10 , L = 5,
D = 0.001 40
15 17
7 Validation du modèle : T en fonction de y, 4> — 0.5,Ra = 10 — 10 ,
L = 5, D = 0.001 41

1 Schéma de fonctionnement de l'algorithme génétique 47

2 Représentation binaire d'un individu, N = 2 49
3 Exemple de décodage d'une chaîne binaire, N = 2 50
4 Exemple de corrections apportées à un individu, N = 5 52
=
5 Sélection des individus par la méthode SUS, -/V,„d 6 et 7Vse;=4 57
6 Croisement de deux individus par la méthode multi-points, m = 3 . . . 58

x
4.7 Schéma du domaine de la recherche locale appliquée au chromosome de
porosités, N = 2, x = 2 60
4.8 Exemple de recherche locale appliquée au chromosome de porosité, iV =
6, x = 1 60
4.9 Exemple de recherche locale appliquée au chromosome de matériaux,
N =b 62
4.10 Exemple de mutation, N = 2 63

5.1 Température maximale dans le stratifié en fonction du nombre de Ray-

leigh et de la longueur,iV = 10, em = ec = 0, D = 0.001 69
5.2 Distribution de porosité dans le stratifié en fonction de x, L — 0.25,
N = 10,e m = ec = 0, D = 0.001 73
5.3 Distribution de porosité dans le stratifié en fonction de x/Leff:Ra, N =
10,e m = ec = 0, D = 0.001 75
5.4 Masse du stratifié en fonction de la longueur et du nombre de Rayleigh,
N = 10, e m = ec = 0, D = 0.001 76
16
5.5 Porosité en fonction de la longueur, Ra = 10 ,e m = ec = 0, D = 0.001 . 77
15
5.6 Structure poreuse optimale sans contrainte, Ra = 10 ,Z = 1.0,iV = 10,
M/Msco = 0, C/Csco = 0, D = 0.001 78
15
5.7 Structures poreuses optimales sous une contrainte de masse, Ra = 10 ,
L = 1.0, iV = 10, Mco/Msœ = [0.75,0.5,0.25,0.1], Cco/Csco = 0, D =
0.001 79
15
5.8 Structures poreuse optimales sous une contrainte de masse, Ra = 10 ,
1=1.0, N= 10, M/Msco = [0.75,0.5, 0.25,0.1], C/C8CO = 0, D = 0.001 82
5.9 Optimisation de la structure poreuse sous une contrainte de masse et de
coût, Ra = 10 15 , L = 1.0, N = 10, D = 0.001 84
5.10 Optimisation de la structure poreuse sous une combinaison de contraintes
de masse et de coût, Ra = 10 15 , L = 1.0, N = 10,M/MSCO = [0.95,0.5,0.2,0.9],
C/Csœ= [0.95,0.5,0.9,0.2], 0 = 0.001 86
5.11 Optimisation de la structure poreuse sous une combinaison de contraintes
de masse et de coût, Ra = 10 15 , L = 1.0, iV = 10, M/Msco = [0.2, 0.3, 0.35,0.35],
C/Csco= [0.1,0.2,0.15,0.3], D = 0.001 87
5.12 Nombre de couches de cuivre dans la structure poreuse optimale selon
la combinaison de contraintes de masse et de côut, Ra = 10 15 , L = 1.0,
N = 10, D = 0.001 88

xi
5.13 Nombre de couches d'aluminium dans la structure poreuse optimale selon
la combinaison de contraintes de masse et de côut, Ra = 1015, L = 1.0,
N = 10,D = 0.001 89
5.14 Nombre de couches d'acier dans la structure poreuse selon la combinaison
de contraintes de masse et de côut, Ra = 1015, L = 1.0, TV = 10,
D = 0.001 90
5.15 Nombre de couches vide dans la structure poreuse selon la combinaison
de contraintes de masse et de côut, Ra = 1015, L = 1.0, N = 10,
D = 0.001 91

XI]
 Liste des symboles

Lettres Latines

c Chaleur latente [j kg"1]

C Coût du stratifié [$]

D Diamètre des pores [m]

g accélération gravitationnelle [9.81 m s~2]

H Hauteur de l'échangeur [m]

k Conductivité thermique [W m'1 K'1]

K Perméabilité [m 2]

L Longueur de l'échangeur [m]

M Masse du stratifié [kg]

N Nombre de couche du stratifié

Nind Nombre d'individus au sein d'une population (AG)

P Champ de pression [Pa]

Pe Nombre de Peclet
//
q flux de chaleur [W m~2]

xiii
Ra Nombre de Rayleigh

Re Nombre de Reynolds

T Température [K]

u, v Composante de la vitesse selon les directions x et y [m s - 1 ]

V Vitesse dans un conduit [m s - 1 ]

x, y Coordonnées cartésiennes [m]

Lettre grecques

a Diffusivité thermique [m2 s~1]

P Coefficient d'expansion thermique [K~l]

ô Épaisseur de la couche limite [m]

e intensité de la pénalité de la fonction F0bjecuf

H Viscosité [kg m2 s" 1 ]

v Viscosité cinématique [m2 s~2]

ijj Fonction de courant

4> Porosité

p Densité [kg m~3]

6 Température [K]

Indices et e x p o s a n t s

c référence au coût

co référence à une valeur contrainte

D référence au conduit

eq référence à une propriété équivalente

xiv
F référence au fluide

i référence à une itération

in référence à l'entrée

m référence à la masse

moy valeur moyenne

num référence à une valeur numérique

out référence à la sortie

s référence au solide

sco référence à une valeur non contrainte

T référence à la Température

theo référence à une valeur théorique

x référence à la direction de l'axe x

y référence à la direction de l'axe y

w référence à une paroi

oo référence au milieu ambient

~ référence à une qantité adimensionnelle

xv
 Chapitre 1

Introduction

1.1 Problématique

Au tournant du 21 e siècle, l'utilisation des circuits intégrés prolifère dans une large
gamme d'appareils de consommation. Notamment, l'effervescence du domaine des té­
lécommunications et des ordinateurs personnels incite les fabricants à concevoir des
composantes électroniques de plus en plus performantes. De plus, d'importants efforts
sont dédiés à la miniaturisation des systèmes. Selon la loi de Moore [1], le nombre de
transistors que contiennent les microprocesseurs sur une puce de silicium double tous
les 18 mois. Depuis les trois dernières décennies, cette prédiction s'est avérée juste.
Cependant, au cours des dernières années, on remarque un certain ralentissement dans
la progression technologique. Les difficultés sont en partie attribuables aux effets quan-
tiques associés à l'infime taille des composants, ainsi qu'aux problèmes de dissipation

I
 2

thermique [2]. En effet, le transfert de chaleur est un aspect fondamental dans le de­
sign des systèmes électroniques modernes [3]. D'ailleurs, le marché mondial du contrôle
thermique est estimé à 5.7 milliards de dollars en 2006 et on juge qu'il grimpera jusqu'à
8.7 milliards de dollars vers 2011 [4]. L'optimisation du contrôle thermique constitue
donc un défi à relever pour l'ingénierie moderne.

Afin de s'assurer de la fiabilité et de la performance d'un appareil électronique, ou

de n'importe quel autre disposif générant de la chaleur, il est primordial d'effectuer une
régulation thermique adéquate. Une des solutions qui permet de maintenir un système
à des températures acceptables est l'emploi d'un refroidisseur. Ce dernier est conçu de
manière à extraire un maximum de chaleur d'une source chaude. Les mécanismes de
transport thermique utilisés pour ce type d'application sont la conduction, la convec-
tion et le rayonnement. Il existe une multitude de types d'échangeur de chaleur, que
l'on retrouve dans une foule d'applications. Dans le domaine de l'électronique, le mo­
dèle le plus courant est de type heat sink. Cet échangeur est généralement constitué
d'une juxtaposition d'ailettes fabriquées d'un matériau conducteur comme le cuivre ou
l'aluminium. Le système se greffe à la source chaude. Par conduction, la chaleur est
diffusée à travers les ailettes, et évacuée en convection à l'aide d'un fluide caloriporteur.
Le fluide employé peut circuler en boucle fermée, comme dans le cas d'un radiateur
d'automobile, ou peut être constamment renouvelé. L'écoulement est soit : induit par
convection naturelle, ou généré par une pompe. Les systèmes qui utilisent la convection
forcée ont la capacité d'extirper une grande quantité d'énergie, mais présentent l'incon-
vient d'être sujets aux bris mécaniques. En contrepartie, les échangeurs qui fonctionnent
par convection naturelle, en raison de leur simplicité inhérente, présentent l'avantage
d'être particulièrement fiables.

En raison de l'accroissement de la charge thermique à laquelle sont soumis les nou­

veaux appareils, les ingénieurs ont tenté d'adapter les heat sinks classiques, constitués
d'ailettes rectangulaires, afin qu'ils puissent dissiper davantage de chaleur. Plusieurs
approches ont été envisagées comme l'optimisation de la forme des ailettes et de la
géométrie du refroidisseur [5]-[6]. Des caloducs (heat pipe) peuvent aussi être ajoutés,
ce qui permet d'accroître la capacité d'extraction thermique [7]. Les matériaux à chan­
gement de phase [8] ainsi que les nanofluides [9] sont aussi envisagés pour maximiser
la performance du refroidisseur tout en réduisant la taille de la structure. L'emploi des
matériaux poreux s'avère également une approche prometteuse [10]. Comparés aux re-
froidisseurs conventionnels, les échangeurs poreux possèdent l'avantage d'offrir un ratio
 3

de transfert de chaleur par unité de volume considérablement élevé. Le fonctionnement

de ce type d'échangeur se compare à celui des systèmes conventionnels. Au sein de la
structure poreuse, un fluide circule à travers le réseau de pores afin d'extraire la chaleur
qui se répartit par conduction dans la partie solide du milieu.

Au cours de cette étude, l'emploi des matériaux poreux dans la fabrication d'un
echangeur de chaleur sera étudié. Plus précisément, il sera question du rôle que joue
l'architecture des pores sur le transfert thermique. Pour ce faire, un stratifié constitué de
la juxtaposition d'une série de différentes couches poreuses sera considéré. La porosité
du milieu, dénotée par 0 dans cet ouvrage, varie entre 0 et 1. Une valeur faible de
la porosité correspond à une architecture qui présente peu de pores. Pour <j) près de
1, le système est très poreux. L'écoulement interne sera mû par convection naturelle.
La porosité ainsi que les matériaux de chaque strate seront optimisés dans le but de
minimiser la température maximale du système, et ce, sous des considérations de masse
et de coût.

1.2 Revue bibliographique

Le transfert thermique en milieu poreux a suscité beaucoup d'intérêt dans le monde

scientifique. Au cours des dernières années, plusieurs chercheurs se sont penchés sur le
rôle que les matériaux poreux pourraient jouer en transfert thermique. Cette section
présente une revue de la littérature sur ce sujet.

Boomsma et al. [10] ont évalué le transfert de chaleur dans un echangeur com­
posé d'une mousse d'un alliage d'aluminium. Le matériau poreux étudié offre un ratio
surface-volume de l'ordre 10 000 m 2 / m 3 , ce qui permet de concevoir des systèmes très
compacts. La porosité du milieu est modifiée en compressant la mousse. Avant com­
pression, la porosité de la structure de type open — cell atteint 95% et possède 40 PPI
(pores per linear inch). Les résultats de l'étude indiquent que la mousse d'aluminium
compressée performe particulièrement bien au plan thermique. En effet, les échangeurs
poreux génèrent des résistances thermiques deux à trois fois moins importantes que les
meilleurs refroidisseurs commerciaux opérant dans les mêmes conditions. Egalement,
Hetsroni et al. [11] arrivent à des conclusions similaires, à la suite d'études expérimen­
tales sur un echangeur poreux fabriqué d'acier inoxydable fritte. Les auteurs démontrent
 4

qu'un refroidisseur constitué d'une structure ayant une porosité de l'ordre de 30% per­
met de dissiper de grandes quantités de chaleur pour un faible volume. Cependant, pour
de faibles porosités, le système requiert une grande puissance de pompage.

Haack et al. [12] ont étudié les propriétés thermiques d'un échanger de chaleur fabri­
qué de mousse métallique. L'expérience démontre que, pour de grandes températures,
le refroidissement est accentué en diminuant la dimension des pores. Cependant, les
mesures indiquent que la diminution de la taille des pores engendre une plus grande
perte de pression. De plus, les auteurs découvrent que pour un écoulement de fluide
calorifique donné, l'augmentation de la densité du milieu accroît le nombre de Nusselt
(le nombre de Nusselt correspond au ratio des échanges thermiques convectifs par rap­
port aux échanges thermiques conductifs). La performance thermique d'une structure
poreuse est donc substantiellement influencée par la géométrie des pores.

L'optimisation de l'architecture interne d'un dissipateur de chaleur est abordée par

Bejan [13]. L'auteur détermine analytiquement, grâce à une approche basée sur l'inter­
section des asymptotes, une configuration idéale de la matière permettant de maximiser
le transfert de chaleur. L'étude suggère des équations analytiques qui caractérisent l'es­
pacement optimal de micro-conduits emplis d'un milieu poreux saturé. Dans un même
ordre d'idées, Muzychka [14] emploie une approche constructale 1 afin de déterminer
l'arrangement interne optimal d'un échangeur de chaleur formé de micro-tubes. Le
principe derrière cette approche est de dimensionner et de positionner des tubes afin
de tirer profit de l'espace de refroidissement, le plus efficacement possible. Les résul­
tats indiquent qu'en accentuant la complexité de la structure, la résistance thermique
globale du système décroît.

Wei et Joshi [16] s'intéressent pour leur part à un dissipateur de chaleur (heat sink)
dédié au refroidissement de composantes électroniques. L'échangeur en question consiste
en un empilement de couches dans lesquelles sont façonnés des micro-conduits (micro —
channel). Pour un flux de chaleur donné, le refroidisseur étudié requiert une puissance
de pompage moindre, car l'empilement présente davantage de surface d'échange qu'une
simple couche. Aussi, l'étude démontre, de manière numérique, que la résistance ther-
1
La théorie constructale de l'optimisation global sous des constaintes locales explique de manière
simple les formes qui émergent dans la nature. C'est l'idée que l'architecture des écoulements provient
d'un principe de maximisation de l'accessibilité des flux, dans le temps, et dans les configurations libres
de changement. Cette théorie fut développée par [15].
 5

mique du système peut être réduite en optimisant la configuration et la longueur des

conduits. Les auteurs découvrent que la juxtaposition de quatre couches peut diminuer
la résistance thermique de 20%. Dans ce même contexte, Jeevan et al. [17] montrent,
grâce à l'utilisation d'un algorithme génétique, que l'optimisation des dimensions des
micro-conduits et qu'un empilement résulte en de meilleures performances de refroidis­
sement. Notamment, leurs résultats suggèrent l'utilisation de canaux plus courts afin
de diminuer la résistance thermique. La géométrie des micro-conduits a été aussi étu­
diée par Foli et al. [18]. Les auteurs déterminent, grâce à une approche analytique puis
d'après un modèle numérique couplé à un algorithme évolutif, les paramètres géomé­
triques des conduits qui favorisent le transfert thermique tout en minimisant la puis­
sance de pompage. L'étude révèle que le rapport largeur/hauteur (aspect ratio) d'une
conduite rectangulaire influence fortement les performances thermiques.

Jiang et al. [19] s'intéressent à deux types de micro-échangeurs, l'un fait de micro­
conduits (MCHE) et l'autre, de micro-canaux poreux (MPHE). Les auteurs étudient
les systèmes autant sous l'aspect numérique qu'expérimentale. Leurs résultats révèlent
que les performances thermiques reliées à l'échangeur MPHE sont meilleures que pour
un système MCHE. Cependant, les pertes de pression associées aux canaux poreux sont
nettement supérieures. Également, les auteurs démontrent que pour les deux types de
refroidisseur étudiés, le transfert de chaleur par unité de volume est extrêmement élevé
comparativement aux échangeurs conventionnels. L'emploi de micro-structure dans la
fabrication d'échangeur de chaleur permet donc de créer des systèmes plus performants
et plus compacts.

En convection naturelle, une relation basée sur des équations de similitude, dévelop­
pée par Bejan [20], caractérise la variation de la température à la paroi verticale d'un
milieu poreux soumis à un flux thermique constant. Les travaux de l'auteur permettent
d'établir un lien analytique entre la distribution de température et les propriétés des
milieux comme la perméabilité et la conductivité équivalente. La structure poreuse étu­
diée est isotrope. L'analyse peut donc s'appliquer à des matériaux tels que des mousses
métalliques.

Déterminer la configuration optimale d'un stratifié poreux soumis à des contraintes

de masse et de coût représente un problème de taille. Pour surmonter cette diffi­
culté, la technique d'optimisation sélectionnée dans ce mémoire se base sur des heuris­
tiques adaptatives, les algorithmes génétiques (AG). Un avantage de cette méthode, est
 (i

qu'elle ne nécessite aucune dérivée ou linéarité de la fonction objectif et des contraintes

pour guider la recherche (voir section 4.2.4). Cette approche, inspirée des mécanismes
qui régissent l'évolution biologique, est employée, notamment dans des domaines tels
que l'économie [21], la gestion [21], la physique [23] et bien d'autres. En ingénie­
rie, on fait fréquemment appel aux AG pour résoudre des problèmes réliés au design
de matériaux composite [24, 25, 26]. D'ailleurs, des opérateurs propres à l'optimisa­
tion de structures laminées furent développés afin d'améliorer les performances de re­
cherche [27, 28]. L'utilisation d'AG dans le design de systèmes thermiques est encore
embryonnaire [29, 30, 31], mais la croissance fulgurante de la puissance de calcul des
ordinateurs permet d'envisager cette approche. En s'inspirant des techniques employés
pour traiter des matériaux composites un AG fut construit afin d'optimiser un échan-
geur thermique. Le détail de la procédure est présenté à la section 4. De plus, en s'ins­
pirant des travaux de [32, 33], une hybridation de l'AG avec une routine de recherche
locale est réalisé pour accroître les performances d'optimisation (voir section4.2.8).

1.3 Objectifs

1.3.1 Objectif principal

Le but de ce mémoire est d'étudier numériquement les phénomènes thermofluides

associés à une juxtaposition de couches poreuses aux propriétés distinctes. L'objectif
est de déterminer la configuration interne optimale du stratifié et de l'écoulement du
fluide refroidissant qui y circule afin de maximiser le transfert thermique. On consi­
dère ici un refroidisseur qui fonctionne par convection naturelle. Le critère choisi pour
quantifier la performance du système est la température maximale atteinte dans la
structure (hot spot). Dans cette analyse, la chaleur transmise au système est dégagée
d'une plaque chaude adjacente au stratifié. L'optimisation sera effectuée à l'aide d'algo­
rithmes génétiques et la résolution des champs de température et de l'écoulement, grâce
à une méthode de volumes finis. L'application pouvant bénéficier le plus directement
des résultats de cette étude est le refroidissement de composantes électroniques.
 7

1.3.2 Objectifs secondaires

Afin d'atteindre le but énoncé à la section précédente, plusieurs objectifs secondaires

doivent être atteints tels :

- Construire un modèle mathématique qui corresponde au présent problème

- Élaborer une solution analytique au problème thermofluide
Elaborer une solution numérique au problème thermofluide grâce aux volumes
finis
- Construire un algorithme génétique qui puisse optimiser la structure sur le plan
de la distribution de la porosité et de l'ordonnancement des matériaux
- Caractériser les structures optimales en fonction des contraintes prescrites
- Juger de l'ensemble du travail effectué.

1.4 Méthodologie et présentation du document

Les spéfications se rattachant au modèle mathématique considéré sont présentées

au chapitre 2 . On y retrouve aussi les équations modélisant le phénomène de convec-
tion naturelle en milieu poreux. Les hypothèses simplificatrices y sont discutées et les
équations sont adimensionnalisées. Les équations différentielles régissant le phénomène
étant très complexe,leur résolution est basée sur un modèle numérique et celui-ci est
présenté au chapitre 3. On y retrouve une présentation des équations discrétisées par
la méthode des volumes finis et la validation du modèle. De plus, afin de s'assurer de
la justesse des résultats numériques obtenus, une étude d'indépendance de maillage
est réalisée, de même qu'une validation en relation avec des résultats présentés dans
la littérature. Le chapitre 3 présente également la résolution analytique des équations
simplifiées en régime couche limite. Le chapitre 4 est entièrement consacré à la méthode
d'optimisation par algorithmes génétiques. Les principaux mécanismes de cette tech­
nique sont traités. Les résultats issus du couplage entre le modèle numérique et l'AG
sont présentés au chapitre 5. Ensuite on retrouve la conclusion et les recommandations.
Au corps de cette étude se greffent une annexe comprennant le script des principales
fonctions Matlab utilisés lors de ce travail (Annexe A).
 Chapitre 2

Modèle mathématique

Cette section du document présente les paramètres géométriques et physiques qui

caractérisent l'échangeur de chaleur formé d'un empilement de couches poreuses. Une
description des composantes du modèle est d'abord donnée (section 2.1). Ensuite, les
équations mathématiques qui régissent le système sont exposées (section 2.2). Pour le
régime couche limite, deux solutions analytiques au problème thermofluide sont pré­
sentées; l'une est issue d'une approximation par ordre de grandeur (section 2.3.1) et
l'autre est obtenue par la méthode intégrale (section 2.3.2). Une section de ce chapitre
porte sur la plage d'opération du modèle (section 2.4). Finalement, on retrouve une
description des matériaux employés lors de l'optimisation (section 2.5).

H
 9

2.1 Présentation de la géométrie du système

Tout d'abord, abordons les principales caractéristiques géométriques et physiques

de l'échangeur de chaleur. La figure 2.1 représente une vue en coupe de ce dernier. La
chaleur est extraite à la paroi située à l'extrémité gauche du schéma. Sur cette frontière,
le flux thermique q" est considéré constant. À la source de chaleur, se greffe un stratifié
constitué d'une juxtaposition de différentes couches de matériaux poreux. L'échangeur
ainsi formé a une hauteur H et une longe totale L. Chaque couche possède la même
épaisseur soit L/N où N représente le nombre total de couches. Cependant, la porosité
et le type de matériau de chaque couche sont variables. Par conséquent, le système
possède 2N variables de design, soit N variables pour les matériaux et autant pour
les porosités. La paroi opposée à la plaque chaude est posée comme étant adiabatique
et imperméable. Par convection naturelle, le fluide refroidissant, initialement à une
température ambiante, dénotée T ^ , pénètre dans l'échangeur par le bas, circule à travers
les pores et ressort à l'extrémité supérieure. La chaleur émanant de la plaque chaude
se diffuse par conduction dans la partie solide du matériau et est ensuite éjectée dans
l'environnement par l'écoulement de fluide caloriporteur. Le point chaud (T max ) se situe
au sommet de la paroi chaude et nous tenterons de minimiser la température qui prévaut
en ce point. L'optimisation permettra d'adapter la structure poreuse afin d'accroître
les échanges thermiques.
 \ sortie

\ J t t t t t, 11 / .
couche 1 couche 2 couche 3 couche 4 ^ -
^ &

iude
a i a 2 a 3 a 4 ^ -s
M K l K2 K3 K4 / .2 H
<h

-e-
(j) 2 <)) 4
ri
Plaq

y °
g

* ^ # ^
y<>

î t t t Vt t. Tt t t
entrée oo

FlG. 2.1 - Schéma de l'échangeur

2.2 Description du modèle mathématique

Dans cette étude, les échanges thermiques reposent sur la dynamique des fluides
en milieu poreux. Un milieu poreux est constitué d'une matrice solide comprenant un
réseau de minuscules cavités reliées entre elles. Ce réseau de pores permet au fluide
de circuler à travers le matériau. À l'échelle des pores, soit au niveau microscopique,
les caractéristiques de l'écoulement sont telles que la vitesse et la pression sont très
irrégulières. Cependant, au plan expérimental, ces quantités d'intérêt se mesurent dans
un volume comprenant plusieurs pores. Sur cette moyenne spatiale macroscopique, ces
quantités d'intérêt varient de manière régulière en considération de l'espace et du temps,
ce qui permet un traitement analytique. Par conséquent, au lieu de résoudre les équa­
tions de conservation de la quantité de mouvement et d'énergie pour la partie fluide
et l'équation de conduction thermique dans la matrice solide, dans chacun des pores,
on procède plutôt à une évaluation moyennée spatialement de l'écoulement et de la
distribution de température. Mentionnons que le modèle mathématique employé dans
ce mémoire requière une certaine simplicité, car lors de l'optimisation par algorithmes
génétiques (voir section 4), le système thermofluide est résolue un grand nombre de fois.

D'abord, on fait l'hypothèse qu'au sein du milieu poreux, il y a équilibre thermody­

namique local. On considère qu'il y a équilibre thermodynamique lorsque localement le
coefficient de transfert thermique entre fluide et le solide est très élevé. Ainsi, en tout
point du stratifié, la température du fluide est la même que celle de la matrice solide :

Tf(x,y)*iTê(x,y) = T(x,y) (2.1)

Il est à noter que bien que la température du solide et fluide aient à la même valeur
localement, à l'échelle du stratifié, un échange thermique se produit entre les deux
médiums. Dans ce document la température symbolisée par la lettre T correspond
donc à la fois à la température du fluide et à la température du solide. Les indices / et
s sont alors éliminés.

La porosité <f> du milieu est définie comme la fraction du volume total occupée par
les pores :

volume des pores

<l> = ; —i— 2.2
volume total

Une couche ayant une porosité égale à 0 signifie qu'elle est composée complètement
d'un matériau solide alors qu'une porosité de 1 correspond à une couche de fluide.
Lorsque l'on traite de convection en milieu poreux, une autre caractéristique entre en
jeu. Il s'agit de la perméabilité. Cette caractéristique du milieu est intrinsèquement liée
à la porosité. La perméabilité, représentée par la variable K, correspond à la capacité
du milieu poreux à laisser circuler le fluide. Cette propriété peut être déterminée empiri­
quement ou à l'aide de corrélations [34]. En général, la perméabilité peut dépendre de la
direction. Dans ce cas, la perméabilité est représentée par un tenseur. Si la perméabilité
est la même dans toutes les directions, il s'agit alors d'un scalaire.

Dans de cette étude, on fait l'hypothèse que les pores se composent d'une multitude
de micro-conduits cylindriques orientés verticalement. Le nombre de capillaires par
unité de volume définit la porosité du milieu. Dans ce cas, la perméabilité de la structure
est anisotrope, soit dirigée uniquement selon l'axe y et nulle selon l'axe x. Le schéma
présenté à la figure 2.2 illustre ce type de structure poreuse.

FlG. 2.2 - Représentation schématique d'un empilement de micro-conduits

Pour un empilement de micro-conduits, la perméabilité peut être estimée d'après la

géométrie du milieu [34] de telle sorte que :

K (0) = -cfiB2 (2.3)

où D correspond au diamètre des conduits. Dans ce type de matériau, la vitesse du

fluide est unidirectionnelle et verticale. Dans ce cas, Kx = 0 et Ky ^ 0.

L'écoulement de fluide au sein du stratifié est considéré comme étant laminaire et en

régime permanent. D'après la loi de Darcy [36], la vitesse d'un tel écoulement dans un
milieu poreux isotrope (Kx = Ky = K) est proportionnelle au différentiel de pression
imposé ainsi qu'à la perméabilité. En présence de la gravité, la loi de Darcy prend la
forme :

K (dP
 L3

où u et v sont les composantes de vitesse dans les directions x et y respectivement, et p

désigne la viscosité du fluide. Dans le cas présent, compte tenu de la structure poreuse
considérée, ces équations se simplifient :

u =0 (2.6)
K (dP \ . .
(27)
» = -7UH
En convection naturelle, l'interaction du champ de gravité avec le gradient de densité
constitue le moteur de la convection (force d'Archimède). L'expansion thermique du
fluide induit donc un mouvement convectif. Pour prendre en considération cet effet,
l'approximation de Boussinesq [35] est considérée dans ce mémoire. Celle-ci stipule
que la variation de densité est proportionnelle à la variation de température. Selon le
premier terme de la série de taylor de cette hypothèse, l'équation d'état caractérisant
la densité est :

p = p0[l-(3(T-T0)} (2.8)

où po est la densité du fluide à une température de référence T0 et j3 est le coefficient

d'expansion thermique du fluide. On souligne que cette approximation demeure valide
lorsque j3AT « 1. Ce qui conduit à :

u = 0 (2.9)

v = --(^ + g(p0-poP(T-T00))\ (2.10)

De l'équation de continuité,V • u = 0, on obtient que dv/dy = 0, car u = 0. Ce qui

implique que la vitesse est seulement fonction de x. Dans ce cas, pour une position x
donnée on a :

^ = -~-P09 + P09P(T-Too) (2.11)

 M

En intégrant l'éq. (2.11) par rapport à y sur la hauteur du stratifié, on obtient l'expres­
sion de la vitesse dans les pores en fonction de x :

L ~Kdy=~L K d y ' ~ m t i + m i 3 L (T-T^dv (2-12)

Sachant que la variation de pression entre l'entrée et la sortie de milieu poreux
correspond à la pression hydrostatique, PH — P0 = —p0gH, il est possible d'éliminer les
termes relatifs à la pression et isoler l'expression de la vitesse moyenne dans les pores.

v(x) = ^^j"(T-T00)dy (2.13)

où T = T(x,y). L'équation (2.13) permet calculer la vitesse dans chacun des pores
grâce à la distribution de température dans le milieux poreux. L'étape suivant consiste
à caractériser l'équation d'énergie dans le système composé d'une partie solide et d'une
partie fluide.

La partie solide du stratifié est constitué d'un matériau conducteur tel que du cuivre
ou l'aluminium. On fait l'hypothèse que la conductivité du fluide refroidissant employé
est beaucoup plus faible que celle du solide. Cette simplification est valide pour un fluide
standard tel de l'air ou de l'eau [37]. Dans ces conditions, on obtient que kf/ks —► 0.
La conductivité thermique équivalente du milieu est alors approximée par la relation
suivante [38] :

keq = ks^-^. (2.14)

où ks est la conductivité thermique du matériau solide et 4> représente à porosité du mi­

lieu (voir éq. 2.2). Cette relation permet de tenir compte de la particularité géométrique
du milieu poreux constitué d'une juxtaposition de micro-conduits disposé verticalement
à travers desquels circule un fluide. L'équation (2.14) caractérise donc la conductivité
équivalente perpendiculairement aux pores. Par conséquent, la formulation employée
au cours de cette étude diffère du modèle de keq habituellement employé pour traiter
des milieux poreux isotropes.
 L5

Considérant l'hypothèse (2.1), la distribution de température au sein du milieu

poreux, composé de la combinaision d'une partie fluide et solide, est gouvernée par
l'équation d'énergie qui s'écrit :

u-VT=—V-(fcegVr) (2.15)

où keq est la conductivité thermique équivalente définie à l'équation (2.14), et pc est la

capacité calorifique volumique du fluide.

En faisant l'hypothèse que la diffusivité selon y est négligeable en raison de l'im­

portance de la convection dans cette direction, l'équation d'énergie simplifiée s'écrit
comme :

vdT = d2aeqT
dy dx2

Les profils de vitesse et de température s'obtiennent grâce à la résolution des équa­

tions couplées éq. (2.16) et éq. (2.13) avec les conditions limites appropriées. La dif-
fusivité équivalente du milieu poreux aeq est définie comme le ratio de la conductivité
équivalente keq (voir éq. 2.14) sur le produit de la densité et de la chaleur spécifique du
fluide :

aeq = —^~ (2.17)

PfCPf

Comme il a été démontré à l'éq. (2.7), la vitesse est uniquement fonction de x.

Cette dernière est donc constante sur toute la longueur d'un pore et adopte un profil
pleinement développé. De plus, il n'y a pas de développement de couche limite de
vitesse à la paroi chaude car chaque conduit est indépendant. Pour cette raison, on
ne peut pas appliquer de condition de non-glissement à cette frontière. Le profil de
vitesse se calcule donc directement à partir du profil de température sans imposer
de conditions limites à l'écoulement à la sortie de l'échangeur. En ce qui concerne
la paroi chaude, le taux de transfert de chaleur demeure constant. La paroi opposée
est considérée comme adiabatique. À l'entrée de l'échangeur, le fluide pénètre à une
 I(i

température T^. Conséquemment, les conditions frontière pour le cas d'un empilement
de micro-conduits sont :

à x =0 ' dx

ax = L ,§£ = 0
ày = 0 ,T = Tœ (2.18)

Afin de permettre une plus grande flexibilité dans l'analyse des équations ther­
mofluides, le problème est traité adimensionnellement. En introduisant les relations
suivantes :

T —T
(2.19)
Hq"/kf
v
* = TT(k
US (2 20)
-
fl {pfcPfH)
1
K{(j)) 1 (D :>,
"M—g = ■&*[* < <2-21>
x = x/H , y = y/H
D = D/H , k =
i

Ra = tlJUl ï- 2.22
kf af fif
âeq = °^ (2.23)
a
f

Le nombre de Rayleigh tel que défini à l'éq. (2.22) ne se base pas sur la perméabilité,
comme il est généralement le cas dans les modèles de convection en milieux poreux
car dans cette étude, la perméabilité est optimisée. La version adimensionnelle des
équations de conservation, éq. (2.13)-éq. (2.16), devient ainsi :

v(x) = K Ra [ f (y) dy (2.24)

 17

Les conditions frontière adimensionnelles décrites aux équations (2.18) s'écrivent

ax=0 à SE = _ l

a x: = L/H SE = n
U
' dx
ày = 0, ,f = 0 (2.26)

2.3 Solution analytique du régime couche limite

Avant d'entreprendre la modélisation numérique (chapitre 3) et l'optimisation de

l'échangeur (chapitre 4), il est souhaitable de connaître l'allure des solutions aux équa­
tions d'énergie (éq. 2.25) et de quantité de mouvement (éq. 2.24). Grâce à des tech­
niques analytiques telles que l'analyse par ordre de grandeur et la méthode intégrale, il
est possible de formuler des relations décrivant la variation de la température à la paroi
chaude du stratifié en régime couche limite. Il faut préciser que le développement théo­
rique qui fait l'objet ce la présente section est le fruit de l'auteur de mémoire. Le modèle
considéré comprend une seule couche ayant une porosité uniforme. Il est à noter que
l'échangeur soumis à l'optimisation par algorithmes génétiques s'avère beaucoup plus
complexe que celui étudié au cours de cette section. De plus, le refroidisseur étudié dans
de ce mémoire n'opère pas nécessairement en régime couche limite, le choix de ce régime
d'opération dans le calcul analytique permet de simplifier grandement les équations de
conservation. Néanmoins, les résultats issus de la démarche analytique pourront servir
de base afin d'estimer la plage d'opération du modèle mathématiques (section 2.4) et
permetteront d'estimer certaines quantités comme l'épaisseur de la couche limite ther­
mique pour un Ra donné. Ainsi, cette analyse indiquera a priori l'ordre de grandeur
des quantités résultant de l'approche numérique.
 18

2.3.1 Approximation par ordre de grandeur du régime

couche limite

Avec l'objectif d'obtenir une approximation au problème thermofluide que repré­

sente l'échangeur poreux, une analyse d'ordre de grandeur est réalisée pour un régime
en couche limite thermique. Cet état correspond à une propagation de chaleur qui s'ef­
fectue dans un milieu semi-infini. L'échangeur possède alors un rapport de forme tel que
la condition limite située à l'opposé de source chaleur (voir figure 2.1) n'influence pas la
distribution de température dans le système. La couche limite thermique représente une
mince région adjacente à la paroi chaude dans laquelle se produit l'essentiel des échanges
thermiques. L'épaisseur de ce domaine, circonscrite par la condition T(x,y)/Toa m 1,
est désigné par ST(y). Pour le cas considéré dans ce mémoire, l'épaisseur ST croît en
fonction de y. Un schéma du développement de la couche limite thermique au sein du
stratifié poreux est présenté à la figure 2.3.

Tout d'abord, débutons l'analyse avec l'équation d'énergie simplifiée éq. (2.16) :

, .dT(x,y) d2T(x,y)

Par ordre de grandeur, on obtient :

l ~ ji < 228 >

y oT
où ST correspond à l'épaisseur de la couche limite thermique. Il est également pos­
sible d'obtenir une expression de la vitesse dans les pores grâce à l'approximation de
Boussinesq éq. (2.8).

K p g 0 ATV
vr *LËJL a (2.29)

où ATy représente l'ordre de grandeur de la différence de température entre l'entrée et

à la sortie du stratifié. On estime que la température du fluide dans la couche limite à la
 19

sortie de l'échangeur est de l'ordre de grandeur de la température à la paroi (Tw (H)).

Le flux thermique à la paroi à y = H s'exprime comme :

AT
f = keq— (2.30)

où AT représente la différence de température entre la paroi à y = H et le mi­

lieu ambiant et ST désigne l'épaisseur de couche limite thermique. En combinant les
éq. (2.28), (2.29) et (2.30), on trouve une expression pour l'épaisseur de la couche
limite thermique en fonction de la position verticale y.

a (2M!
' w~^r
Grâce aux équations (2.31) et (2.30), il est possible d'estimer la température en
régime couche limite sur la paroi chaude en fonction de la position verticale. Ainsi, AT
devient :

Fto-Too)-!?,01"™) (2.32)
\keq gppKJ

Sous la forme adimensionnelle l'expression (2.32) s'écrit :

fw (y) ~ ( i-*r-\ (2.33)

Ainsi, par une analyse d'ordre de grandeur, on découvre que la température à la paroi
chaude devrait augmenter selon Tw (y) oc y1/3. En raison du développement vertical de
la couche limite thermique, il est possible de déduire que la température maximale
du système se situe à la paroi chaude dans la partie supérieure de l'échangeur comme
l'illustre la figure 2.3. De plus, l'analyse par ordre de grandeur prédit que la température
maximale dans le milieu poreux sera proportionnelle à :
 20

1/3

■L max ~ (2.34)
KRa keq K t

T 4 A u t.
T A sortie
V
I t M4

hrff t M *
FiG. 2.3 - Développement de la couche limite thermique et température maximale dans
le stratifié

2.3.2 Approche par la méthode intégrale du régime couche

limite

Une approche basée sur la méthode intégrale est utilisée afin de raffiner les résultats
analytiques obtenus grâce à une analyse d'ordre de grandeur. Cet exercice a aussi pour
objectif de caractériser la variation de température à la paroi chaude en régime couche
limite.

D'après la loi de Fourier le flux thermique à la paroi gauche se définit comme :

q" = —keq^~. En négligeant la conduction dans le sens de l'écoulement, l'équation
d'énergie (2.16) prend la forme suivante :
 21

f ^ = «'■eq
f (2.35)
^/y dx2

Pour une position y donnée, on intégre selon x de 0 à X, où X est juste un peu à

l'extérieur de la couche limite.

d(T{x,y)-T« d(T(x,y) -Toc)

(Y (\ eq
Jo oy eq dx dx
x=X J x=o
(2.36)

L'expression précédente se simplifie grâce aux conditions limites du système. À la

frontière chaude, le flux thermique est constant, — keq—g~^ = q", et à la frontière
opposée, posée adiabatique, le flux thermique est nul. Dans ce contexte, l'équation 2.36
s'écrit :

= d fx
^~ v (T(x,y)-T00)dx (2.37)
pcp oy Jo

Avec la méthode intégrale, il est nécessaire de poser un profil de température réaliste

afin de poursuivre l'analyse. Au sein d'une couche limite thermique (x < ô?) formée
par la convection naturelle sur une plaque disposée verticalement, la température est
maximale près de la paroi chaude et diminue progressivement en s'éloignant de la source
de chaleur. La variation de vitesse évolue de manière similaire à la distribution de
température. On pose alors les profils de vitesses et de température suivants :

v (x) VI- x
r
x
T(x>y)s-Toc (TUîrt-rooHi (2.38)

La constante V représente la vitesse de l'écoulement très près de la paroi chaude

et Tw correspond à la température en fonction de y à la paroi. Les équations (2.38)
sont valables à l'intérieur de la couche limite thermique soit pour x < ÔT- En évaluant
 22

(2.37) d'après les expressions (2.38), on obtient une relation permettant de connaître
l'épaisseur de la couche limite thermique 8T en fonction de la position verticale y :

'3r> \l> 2
! -%h) (2-39)

Selon cette analyse, l'épaisseur de la couche limite croît selon y1/2. La vitesse ver­
ticale du fluide dans le milieu poreux est calculée grâce à l'équation (2.13). À la paroi
chaude, l'expression de la constante V est :

v = 9P PK f
/ (Tw - Too) dy
Vo

En introduisant la relation Tw — T ^ = q"6T/keq, la vitesse V s'écrit alors comme ;

y gf3pKq" rïîr"(3a
/ 3 ne \ ll'2
fiH keq Jo \ V

Sous la forme adimensionnelle, la vitesse à proximité de la paroi chaude est évaluée

d'après le résultat de l'intégrale précédente. Le résultat conduit à :

V = — ^ = 1.006 (Ra K ~keq)2/3 (2.41)

cxf/H ^ '
où V représente la forme adimensionnelle de la vitesse de l'écoulement près de la paroi
chaude définie à l'éq. 2.38. Il est alors possible de déduire la répartition de température
à la paroi chaude du stratifié.

* W [H) J- OO - , n TJ '■ TT ~- TJ \ T/ V
k q„H H H\ V

Sous la forme adimensionnelle, ce résultat s'écrit :

 2:5

/ i v 1/3
T 1M f m <242)
-®- {l*lZlt) '
Suite à cette analyse par la méthode intégrale, on déduit que la température à la
paroi chaude varie proportionnellement à y1/2. La température maximale atteinte au
sein du stratifié, située ky = H sur la paroi chaude, diminue avec le nombre de Rayleigh
selon :

L/3
T
- -L65 U Ï ^ J (2 43)
-
De plus, l'épaisseur de la couche limite à y = 1 correspond à :

/ ~e \1/3
*r = 1.65 \-^f-^\ (2.44)
\RakeqK I

Ces dernières relations sont similaires aux résultats issus de l'approche par ordre de
grandeur. En comparant les équations éq. (2.33) et éq. (2.42), on remarque que l'aug­
mentation de température à la frontière chaude, dans les deux cas, est proportionnelle
à Ra~1^ k~qlz K~1/3. Cependant, l'expression (2.42) suggère que la température varie
selon y 1 / 2 , alors que l'équation (2.33) propose une variation plutôt selon y1/3. Cette
différence entre les exposants provient de l'estimation de la vitesse dans la couche li­
mite. Lors de l'analyse d'ordre de grandeur, on estime que la vitesse demeure constante
selon x alors que pour le cas de la méthode intégrale on propose plutôt un profil de vi­
tesse décroissant (éq.2.38). Les résultats du modèle numérique, présentés ultérieurement
(chapitre 3), démontreront que l'approximation issue de la méthode intégrale s'avère
être plus représentative de la distribution de température à la paroi. Néanmoins, les
deux techniques employées concordent quant à l'estimation de la valeur et de la position
de la température maximale dans l'échangeur.
 24

2.3.3 Porosité optimale en régime couche limite

Grâce aux expressions analytiques développées dans les sections précédentes, il est
possible d'estimer le profil de température à la paroi chaude en fonction des paramètres
du fluide refroidissant, de l'intensité de la sollicitation thermique et des propriétés des
couches du stratifié. D'après les résultats de l'approche intégrale, la température maxi­
male dans l'échangeur est (éq. (2.43)) :

( i Y 1/3
T m a x = 1.65 =——
\Ra keq K)

L'objectif principal de cette étude est de déterminer l'arrangement optimal des

couches poreuses afin de minimiser la température T max . En introduisant les expres­
sions de la conductivité équivalente (éq.( 2.14)) et de la perméabilité (éq.( 2.3)) dans
l'équation (2.43), on obtient :

Ainsi, pour un fluide donné et un échangeur doté d'une seule couche poreuse aux
propriétés constantes, on peut exprimer la température maximale en fonction de la
porosité. La fonction (2.45) évaluée entre 0 et 1 atteint un minimum situé à 4>opt = 0.414.
Cette valeur correspond à une porosité optimale, en régime couche limite pour une seule
couche.

2.4 Plage d'opération du modèle

Selon les hypothèses posées lors de la modélisation mathématique (section 2.2) on

stipule que la diffusité thermique selon y peut être négligeable dans l'équation d'énergie
éq. (2.16). Cela requiert que dans cette direction, les échanges reliés à la convection
dominent par rapport aux effets conductifs. Ce rapport peut s'exprimer par le biais du
nombre de Peclet. Ainsi, pour respecter l'hypothèse il faut que :
 25

V D
Pe = >> 1 (2.46)

De plus, l'équation de Darcy tel que présenté à l'éq. (2.5) n'est valable que dans
le cas d'un écoulement laminaire. Pour un écoulement en conduite, la transition de
laminaire à turbulent s'effectue autour d'un nombre de Reynolds de 2300 [39]. Ainsi,
l'écoulement dans un pore doit respecter :

V D
ReD = < 2300 (2.47)

où V est la vitesse moyenne maximale dans un pore situé près de la plaque chaude et D
le diamètre du pore. La vitesse calculée selon la loi de Darcy correspond à une vitesse
moyennée. La vitesse réelle dans les pores est :

V=~ (2.48)

En substituant l'expression adimensionnelle de la vitesse et du diamètre dans l'évalu­

ation du nombre de Reynold, on obtient la relation suivante :

V D
ReD = — — < 2300 (2.49)
6 Pr

En outre, la vitesse de l'écoulement à proximité de la paroi chaude peut être évaluée

grâce à l'expression développée par la méthode intégrale éq. (2.41). La relation éq. (2.49)
prend alors la forme de :

1.006 RaWRWkWÙ , N
ReD = %— - < 2300 2.50
ç Pr
où le nombre de Prandtl correspond à Pr = v/a. De cette dernière expression, éq. (2.50),
on conclut que le nombre de Reynolds dans les micro-conduits est fonction des pa­
ramètres géométriques du stratifié par l'entremise des variables <j> et D, du type de
matériau qui constitue l'échangeur avec keg, du nombre de Rayleigh et du nombre de
 20

Prandtl. Dans le cas où le fluide refroidissant est de l'air ou de l'eau, le nombre de

Prandlt demeure autour de Pr ~ 1. De plus, pour un echangeur constitué de métaux
conducteurs (tab. 2.1) et ayant une dimension de pore égale à D = 0.001, la plage de
nombre de Rayleigh qui respecte les conditions éq. (2.46) et éq. (2.47) est :

10 13 < Ra < 10 19 (2.51)

Les simulations réalisées en dehors de cette plage de Ra nécessiteront des modèles

mathématiques et numériques beaucoup plus complexes afin de tenir compte des parti­
cularités de l'écoulement et du transfert de chaleur pour de tels nombres de Rayleigh.

2.5 Propriétés des matériaux employés

L'objectif de cette étude est de déterminer la configuration optimale du stratifié

poreux de manière à minimiser T m a x . Selon certaines conditions d'application, il est
souhaitable de contraindre l'echangeur optimal à respecter une masse et/ou un coût
spécifique. Lors de cette analyse, seulement des métaux sont considérés pour constituer
la matrice solide du milieu poreux. Les propriétés de ceux-ci sont énumérées au tableau
(2.1). Les valeurs de conductivité et de densité présentées sont normalisées avec les
propriétés de l'air évalué à 300 K [37]. Les coûts des matériaux sont normalisés par
rapport au prix de l'acier [40]l.

TAB. 2.1 - Propriétés des ma tériaux disponibles lors de l'optimisation

Matériaux Densité, p Conductivité, k diffusivité, à Prix, c
Alluminium (al) 2327.5 9011.4 5.20 12.8
Cuivre (eu) 7691.6 15247.2 4.32 18.9
Acier (st) 6776.3 3049.4 1.03 1
Vide (vi) 0.0 0.0 0.0 0.0

Considérant que la masse de l'air est beaucoup faible que la masse des métaux
énumérés au tableau précédent, celle-ci est négligée et la masse totale du stratifié est
^ e s propriétiés thermiques de l'acier correspondent à celles du fer pur qui représente la limite
supérieure en terme de conductivité thermiques
 27

supposée égale à la somme des masses individuelles de chacune des couches. Ainsi, la
masse normalisée se calcule de la manière suivante :

N
M
M'
= ^r-m=EpAl-<Pj) ^J (2-52)
Pîn j=i

où Aiij, pj et (f>j sont respectivement l'épaisseur, la densité et la porosité de la j e couche

et M représente la masse de l'échangeur. Similairement, le prix de l'échangeur peut être
calculé en se référant au prix des matériaux qui composent chaque couche. Le coût total
d'un échangeur est calculé d'après :

C = Y,m(1-(f>j)A£j (2-53)

e
où ëj représente le coût normalisé de la j couche.

2.6 Conclusion

L'optimisation d'un échangeur constitué d'une série de couches poreuses juxtaposée

dans lesquelles le fluide refroidissant circule par convection naturelle (figure 2.1) vise
à minimiser la température maximale ' atteinte à la plaque chaude. Pour ce faire, un
modèle mathématique est développé. Le type d'architecture poreuse abordé est une
juxtaposition de micro-conduits. La température maximale du système est fonction
des caractéristiques des couches du stratifié (k, (f>, Ra, D, L, N). On estime la valeur de
Tmax, grâce à une approximation par ordre de grandeur (éq. 2.33) et par la méthode
intégrale (éq. 2.42). Ces résultats analytiques sont ensuite employés pour déterminer
la plage d'opération du modèle. Toutefois, les équations d'énergie et de quantité de
mouvement (éq. (2.24) et éq. (2.25)) présentées dans ce chapitre devront être résolues
numériquement. Le chapitre suivant traite du modèle numérique employé pour réaliser
cette tâche.
 Chapitre 3

Modèle numérique

La méthode des volumes finis se révèle être une approche très efficace dans la réso­
lution de systèmes sujets à des lois de conservation. Elle est d'ailleurs largement utilisée
en mécanique des fluides [41]. Cette méthode sera donc employée dans cette étude pour
résoudre le problème thermofluide présenté au chapitre précédent. Au cours du pré­
sent chapitre, les principaux éléments du modèle numérique seront abordés tel que : le
maillage utilisé (section 3.1 ) et la discrétisation du domaine (section 3.2 ). Ensuite,
l'indépendance de maillage et la validation des résultats seront traités (section 3.3 et
section 3.4).

'AH
 2!)

3.1 Maillage utilisé

Le domaine de calcul est fractionné uniformément en cellules rectangulaires sur les­

quelles les équations d'énergie, éq. (2.25), et de conservation de mouvement, éq. (2.24),
sont résolues. La figure 3.1 illustre les principales caractéristiques du domaine et le
maillage utilisé. On dispose Ny cellules dans la direction y. La même densité de cellule
est utilisée dans la direction x. Le nombre total de volumes de contrôle correspond alors
à NyL où L est la largeur de l'échangeur comme l'indique la figure 3.1. Le nombre de
couches physiques N qui composent le stratifié est un paramètre de l'optimisation. Ces
couches physiques sont fractionnées en sous-couches de cellules numériques ayant des
propriétés identiques. À la figure 3.1, les cellules d'une même couleur identique appar­
tiennent à la même strate physique. Dans cet exemple, le stratifié se compose donc de
quatre couches. Afin de respecter la densité de cellules Ny, le nombre de sous-couches
par couche physique équivaut à Ny L/N. En raison du type de discrétistion employé,
l'épaisseur des mailles aux frontières de l'échangeur possède une superficie égale à la
moitié des cellules du domaine central.

mtMmmHM
-xi
Kl

eu
/ s
7-1
o' y -
ai
s P
s "-S

y n

mmmmtm
f 4-
3 4
.

couche
physique
V. T = 0

FlG. 3.1 - Exemple de maillage : L = 1, N = 4, Ny = 16

 30

3.2 Discrétisation par volumes finis

Cette méthode de calcul numérique repose sur des bilans de masse, de quantité de
mouvement et d'énergie sur des volumes de contrôle qui sont définis par le maillage.
La méthode des volumes finis permet de remplacer des équations différentielles par
un système d'équations algébriques. Chaque noeud du maillage est situé au centre
d'un volume de contrôle élémentaire. Sur les frontières de ce volume, les équations
différentielles partielles sont converties en intégrales de contour. Les gradients s'évaluent
alors comme des flux. Le flux qui entre dans un volume est identique à celui qui sort
du volume adjacent, il y a donc conservation. Pour un système bidimensionnel, les
volumes se réduisent à des rectangles de contrôle. Comme l'illustre la figure 3.2, dans de
cette étude, un maillage structuré est employé. Les équations qui régissent le problème
thermofluide sont évaluées autour des points situés au centre des cellules. Une notation
inspirée des points cardinaux sert à identifier les cellules entourant le volume central.
Ainsi, les indices n, s, e, w désignent respectivement les positions situées nord, au sud,
à l'est et à l'ouest du volume central indiqué par la lettre p. À la figure 3.2, les indices
inscrits en lettre majuscule correspondent aux centres des cellules voisines, alors que
les indices minuscules représentent les différentes frontières de la cellule centrale.

U
N

M * * * ** K
n
y
.ZZ-Y,
x- u
w I LE
"- - 7 S •
C> '.115
.. .
EZ2
..-',/'
y
y
y
T
w#
W
K

w
e L
E

:::: / K
::__ y,
:: : /, s
y U
• S
kM M M Ts

FlG. 3.2 Discrétisation par volumes finis

Selon cette méthode de discrétisation, l'équation d'énergie (2.16) évaluée sur les
mailles situées à l'intérieur des conditions limites devient :
 31

La vitesse au centre des surfaces de contrôle est calculée d'après un schéma amont
du premier ordre. Ainsi, l'amplitude de la vitesse à la limite du volume de contrôle d'une
cellule s'estime grâce à l'intensité de la vitesse de la maille située en amont. En considé­
rant ce schéma, de même que la formulation de la vitesse démontrée à l'équation (2.13),
le terme relatif à la convection dans l'équation d'énergie devient :

te i-n Q

L h dy ^ (X) T (X'V)] dy dX = A'X [V (X)" Tn ~ V (x)s Ts] (3 2)

'
Dans un milieu poreux constitué de micros conduits, la vitesse est constante pour
une position x donnée. Ce qui mène à la simplification suivante : v (x)n = v (x)s = v (x).
Le terme convectif de l'équation d'énergie peuvent alors s'exprimer directement en
fonction des températures et des vitesses calculées aux points P et S, comme suit :

f fnv (x) dT( V

f' ^dy dx = Axv (x) (TF - Ts) (3.3)
Jw Js ' dy

La même procédure est appliquée aux termes de diffusion. La conductivité équiva­

lente, Keq, s'évalue à la frontière de la surface de contrôle selon une moyenne harmonique,
telle qu'indiquée à la figure 3.2.

*d ( dT(x,y)\^^ _ A f fdT(x,y)\ fdT(x,y)\\

dxdy = A y
dx\^-dx— [ ^ [ - ^ x ~ - ^ M -dx
(3.4)

s:si (- *¥) ** - » ( - ^ - - ^ ) ™
Une fois que la discrétisation des équations différentielles est terminée, le système
 32

peut alors s'exprimer sous une forme algébrique telle que :

aP TP + aN TN + as Ts + aE TE + aw Tw = r

Les coefficients a et le terme source r s'expriment comme :

aN 0 (3.6)
as = Axv(x)
Ay Keqw
aw Ax
Ay KeQe
aE
Ax
r -asTs (3.7)

La température en chacun des points du domaine est calculée grâce à la résolution

du système AT = r, où la matrice A, formée des coefficients a, est tridiagonale.
ap aE 0 0 0... "Ti" n
a\y ap aE 0 0... T2 T2
= (3.i
0 a\y ap aE 0...

Le processus de résolution est itératif. Une répartition de température est d'abord

calculée d'après un profil de vitesse donné. Les équations de quantités de mouvement
et d'énergie étant couplées, on obtient un nouveau profil de vitesse correspondant à la
répartition de température calculée. Ensuite, ce nouveau profil de vitesse est employé
dans une deuxième résolution de l'équation d'énergie. Ces opérations se répètent jusqu'à
ce que la convergence soit atteinte. Lorsque les normes des changements relatifs du profil
de température et de vitesse entre une itération et la suivante sont inférieures à 0.1%,
on considère que la convergence est obtenue.

1 ri ( fijxM-fi-^xâ) dx dy < 10"

Jo Jo i(x,v)
ri fvi(x)-Vi-i(x)\ i~
ax < 10" (3.9)
h { St-iC*) S

où (i) est l'itération en cours.

 33

3.2.1 Exemple de résultats numériques

Une fois que la discrétisation par volumes finis est complétée, il est alors possible
de calculer la distribution de température et de vitesse donnée. Un exemple de résul­
tats obtenus grâce au modèle numérique décrit dans le présent chapitre est illustré à
la figure 3.3. Il s'agit d'une simulation réalisée selon les paramètres recommandés par
l'indépendance de maillage, traitée à la prochaine section 3.3. L'image de gauche in­
dique la distribution de température au sein de l'échangeur. On remarque nettement
le développement de la couche limite thermique à la paroi. Ce phénomène se manifeste
par mince région colorée adjacente à la source de la chaleur. La couleur variant du bleu
au rouge souligne le changement de température qui se produit dans le milieu poreux.
La température maximale se retrouve au sommet de la paroi, comme il a été démontré
à la section 2.3 et illustré à figure 2.3. La portion de droite de l'illustration présente
la distribution de vitesse dans les pores du refroidisseur. La vitesse maximale s'observe
près de la source de chaleur. Le fluide initialement à une température froide, se réchauffe
au contact de la structure solide chaude, ce qui induit un écoulement vers le haut dans
les micro-conduits. En s'éloignant de la paroi, les effets convectifs s'amenuisent rapide­
ment. On constate alors que pour la longueur d'échangeur considérée L — 5, une large
partie du milieu poreux ne contribue pas aux échanges thermiques.

FlG. 3.3 - Distribution de température et profil de vitesse dans un échangeur poreux,

0 = 0.5, Ra = 10 15 , 1 = 5, D = 0.001

Un deuxième exemple des résultats issus de simulations thermofluides est illustré à la

figure 3.4. Cette fois, une courte longueur d'échangeur est choisie (L = 0.5). On note que
la chaleur se répartit sur l'ensemble du milieu poreux. En effet, le système n'évolue pas
 M

en régime couche limite. D'ailleurs, la vitesse au sein de la structure demeure de forte

intensité dans la totalité des pores. Finalement, pour cette longueur, la température
maximale dans le refroidissement est 17% plus élevée que dans le cas opérant en couche
limite thermique.

0.2 0.3 0.4 0.5

FlG. 3.4 - Distribution de température et profil vitesse dans un échangeur poreux,

0 = 0.5, Ra = 10 15 , L = 5, D = 0.001

3.3 Indépendance de maillage

Une fois le modèle numérique complété, une analyse d'indépendance de maillage

s'impose afin de s'assurer de la validité des résultats obtenus. L'étude d'indépendance
s'effectue sur un échangeur constitué d'une seule couche ayant une porosité donnée.
Comme il a été décrit à la section 3.1, le nombre de noeuds par unité de longeur Ny est
uniforme sur l'ensemble du domaine. Ainsi, une première simulation numérique est réa­
lisée avec un nombre de cellules donné. Du profil de température calculé, la température
maximale est extraite. Ensuite, pour les mêmes conditions d'opération, un deuxième
calcul est réalisé avec une densité de nœuds doublée. Une nouvelle température maxi­
male est alors calculée. Le maillage est ainsi raffiné jusqu'à ce que la variation relative
de la température maximale au cours des raffinements successifs soit inférieure à 1 %.
Ce critère indique donc l'indépendance de maillage.

maxj max,j+l
<i% (3.10)
:/;max,j-
 35

où (j + 1) est le maillage en cours. Le tableau 3.1 présente les résultats d'analyse

d'indépendance de maillage obtenu pour un stratifié fait d'une seule couche d'aluminium
ayant une porosité de 0.7.

T A B . 3.1 - Étude d'indépendance de maillage pour un stratifé fait d'aluminium, L — 1,

0 = 0.7, N= 1, D = 0.001
Ra = 1014 Ra = 1018
Ny f
J
■*max,i -'max.z + l
J- max,t+l
f
x
-* max,i -* max,i f 1

max max
15 1.4568E-3 - 4.1717E-5 -
30 1.4546E-3 0.15 % 4.4871E-5 7.56 %
60 - - 4.5502E-5 1.41 %
120 - - 4.5876E-5 0.82 %

La même procédure est répétée pour les différents matériaux soumis à l'étude. Après
l'analyse des résultats numériques, on constate que la densité du maillage n'est pas
influencée par le type de matériau constituant le refroidisseur poreux. Cependant, il faut
préciser que cette conclusion s'applique seulement aux métaux conducteurs énumérés
au tableau 2.1. De plus, cet exercice révèle que pour un matériau donné, la densité de
noeuds requise est fonction à la fois du nombre de Rayleigh et de la porosité <fi. Pour des
porosités élevées, l'étude d'indépendance recommande un plus grand nombre de cellules
que pour de faibles porosités. Une comparaison des résultats montrés aux tableaux 3.1
et 3.2 démontre bien la nécessité d'affiner le maillage avec l'augmentation de la porosité.
Ces tableaux présentent la variation de température maximale en fonction de la taille
du maillage pour différentes valeurs de 0 d'un système donné. On constate aussi que la
densité de noeuds Ny requise augmente en fonction du nombre de Rayleigh. En d'autres
termes, on retrouve le maillage le plus fin pour un Ra et un (f) élevés. Cette relation est
présentée au tableau 3.3, où l'on retrouve la variation du Ny requis en fonction du Ra
et de (j). Egalement, des stratifiés composés de divers matériaux aux porosités variées,
semblables à ceux présentés au chapitre 5, ont été étudiés. L'examen indique que la
densité de noeuds requise pour atteindre l'indépendance de maillage dans le cas d'un
échangeur fait de plusieurs couches différentes est moindre que celle exigée pour une
structure dotée d'une porosité très élevée (0=0.95) et ce pour la gamme de nombre de
Rayleigh considérée.

À la lumière de ces résultats, on note que le nombre de Rayleigh influence grande­

ment la discrétisation du domaine, et par le fait même, le temps de calcul. Afin d'opti-
 36

T A B . 3.2 - Étude d'indépendance de maillage pour un stratifé fait d'aluminium, L

0.5, 0 = 0.95, N = !,£> = 0.001
Ra = 10 14 Ra = 10 18

Ny f1 -*max,i ^max,i-fl
-* max,i + l
f± -*max,i Jmax,i|l
■*■ m a x , î + l
- max max
15 l,7430E-3 - 4,657E-5 -
30 l,7459E-3 1.66 6,148E-5 32,02 %
60 l,7462E-3 0.017 6,881E-5 11,91 %
120 - - 7,598E-5 10,43 %
240 - - 7,825E-5 2,98 %
480 - - 7,884E-5 0,75 %

T A B . 3.3 - Ny requis pour l'indépendance de maillage en fonction du Ra et de <f> pour

un stratifé fait d'aluminium, L = 0.5, iV = 1, D = 0.001
0 Ra = 10
14
lia - 10" Ra = 1016 Ra = 1017 Ra = 1018
0.10 30 30 30 30 30
0.20 30 30 30 30 30
0.30 30 30 30 30 60
0.40 30 30 30 30 60
0.50 30 30 30 30 60
0.60 30 30 30 30 60
0.70 30 30 30 60 120
0.80 30 30 30 60 240
0.85 30 30 60 120 240
0.90 30 30 60 240 480
0.95 30 60 120 240 480
 37

miser l'algorithme, le nombre de noeuds employés pour résoudre le système s'ajuste en

fonction du nombre de Rayleigh. Lors du processus d'optimisation, des porosités variant
entre 0 et 1 sont considérés, et le maillage le plus fin sera employé indépendamment de
la porosité. Par conséquent, les valeurs de Ny sélectionnées par la suite correspondent
à la dernière ligne du tableau 3.3.

3.4 Validation du modèle numérique

La validation du modèle s'effectue en comparant les résultats numériques issus de

la résolution des équations d'énergie et de conservation de la quantité de mouvement,
présentés à la section 3.2, avec des résultats provenant du logiciel commercial de modéli­
sation thermofluide Fluent [42]. De plus, le modèle analytique développé à la section 2.3
est employé pour consolider la validation du code CFD. Lors de l'analyse, l'échangeur
considéré est constitué d'une seule couche ayant une porosité uniforme. La validation
est réalisée en régime d'écoulement de type couche limite. Pour ce faire, la longueur L
de l'échangeur, illustrée à la figure 2.f, est suffisamment importante pour que la couche
limite thermique ainsi que la couche limite de vitesse ne soient pas influencées par la
paroi adiabatique à la sortie de l'échangeur soit à y = f. Numériquement, ce critère est
déterminé selon :

T(x = L, y= l)
-J: f < 1% (3.11)
T(x = 0, y=l)

où T (x = L, y — 1) est la température à la frontière adiabatique et T (x = 0, y = 1),

la température maximale. Lorsque ce rapport est inférieur à 1%, on fait l'hypothèse
que l'écoulement à \x = L, y = 1 ne ressent pas la présence de la plaque chaude. Au
tableau 3.4, on retrouve les valeurs de L requises pour répondre au critère (3.11) dans le
cas d'un milieu ayant des porosités de 4> — 0.5 et <>
/ = 0.9. Conformément au résultat de
l'approximation obtenue par la méthode intégrale (éq. 2.44 ), l'épaisseur de la couche
limite thermique décroît en fonction du nombre de Rayleigh.

Afin que les profils de température et de vitesse calculés sous Fluent reflètent les don­
nées du modèle numérique de ce mémoire, les propriétés du fluide refroidissant ainsi que
 38

T A B . 3.4 - Longeur requis pour un développement couche limite en fonction du Ra,

(f> = 0.5 et 4> = 0.9 , N = 1,D = 0.001
Ra 1014 1015 1016 1017 1018
0 = 0.5 L 4.21 1.95 0.90 0.42 0.19
0 = 0.9 L 1.01 0.47 0.21 0.10 0.047

du matériau fournies au logiciel sont adimensionnalisées. De plus, conformément à l'hy­

pothèse présentée à l'éq.(2.16), la conductivité dans le sens de l'écoulement est négligée.
Une perméabilité unidirectionnelle correspondant à un empilement de micros conduits
est aussi imposée au modèle. Cependant, le calcul de la conductivité équivalente effec­
tué par Fluent diffère de celui employé dans cette étude ( éq. (2.14) ). L'expression de
keq selon Fluent est :

keq = ks(l-(f>) + kf(p (3.12)

où ks représente la conductivité du solide et kf la conductivité du fluide. Par conséquent,

pour garantir la similitude entre les modèles, l'équation (3.12) est utilisée temporaire­
ment lors des simulations numériques destinées à être validées. Les mêmes conditions
limites qu'énoncées à la section 2.2 sont imposées aux parois chaudes et adiabatiques.
Par contre, sous Fluent, il est difficile de fixer des conditions limites identiques à celles
utilisées dans ce mémoire, à l'entrée et à la sortie de l'échangeur. En effet, étant donné
que le milieu poreux est constitué d'un empilement de micros conduits, la vitesse à l'en­
trée varie selon la position en x par rapport à la plaque chaude, comme il est expliqué à
la section (2.2). L'amplitude de la vitesse du fluide dans un conduit est liée à l'intensité
des échanges thermiques qui s'effectuent dans le tube selon la direction y. Conséquem-
ment, à l'entrée ainsi qu'à la sortie, l'écoulement est dirigé verticalement et distribué
non uniformément selon x. Ces particularités ne peuvent pas être facilement imposées
dans Fluent. Pour palier à ce problème, deux régions virtuelles supplémentaires sont
ajoutées aux extrémités du milieu poreux étudié. Une première région poreuse se greffe
à la sortie de l'échangeur. Celle-ci possède une grande perméabilité selon y et une très
faible perméabilité en x. Cet artifice permet de forcer l'écoulement à demeurer vertical
à la sortie sans imposer de conditions supplémentaires qui pourraient affecter les résul­
tats. De plus, la grande perméabilité verticale réduit les pertes par friction associées
à l'ajout de matière. Une deuxième section constituée uniquement de fluide (0 = 1)
 39

s'annexe à l'entrée créant ainsi une zone tampon entre le stratifié et la condition li­
mite imposée dans Fluent. Cette dernière modification a pour objectif de permettre au
fluide de pénétrer librement dans le milieu poreux sans être influencé par les conditions
imposées à l'entrée du système. Grâce à l'addition de ces zones, il est possible de fixer
des conditions de pressure inlet et pressure outlet à l'entrée et à la sortie du modèle
Fluent. Le différentiel de pression imposé est nul de sorte que le moteur de l'écoulement
soit uniquement la convection naturelle. La figure 3.5 illustre le modèle formé du milieu
poreux central et des deux régions supplémentaires. Le maillage du système, réalisé
avec le logiciel Gambit [43], respecte les caractéristiques décrites à la section 3.1 quant
aux nombre de noeuds linéique Ny et à la forme des cellules.

Pressure outlet - 0

I- *i
4 4 4 k\i 44 4 4 4

zone : sortie
K = 0, Ky= »
6 = 0.95

milieu poreux
K = 0, K=l/32d>D2
6 = 0.5
l H

/ i ■-

zone : entrée
K = oo, K = oo
6=1

JA
L i A A
V entrée
7^
' oo
AA A A
Pressure inlet ■* 0

FlG. 3.5 - Exemple du maillage utilisé avec le logiciel Fluent : L = 1, N = 1, Ny = 16,

D = 0.001

Pour différents paramètres d'opération, mais toujours en régime en couche limite,

on compare le profil de température à la paroi chaude obtenue avec le logiciel Fluent à
celui issu du modèle numérique présenté à la section 3.2. L'analyse d'ordre de grandeur
et l'approche par la méthode intégrale révèlent que la température à la paroi évolue de
comme une puissance de y. Ainsi, l'équation (2.42) indique que :
 40

Tw (y) = ay;. b (3.13)

où a = Ra~]/3 k~q^ K'1^ et b = 1/2. Bien que cette relation soit approximative,
les résultats numériques doivent néanmoins refléter une certaine similitude. La figure
(3.6) présente les résultats obtenus pour le profil de température à la paroi chaude d'un
échangeur poreux en aluminium doté d'une porosité de 0 = 0.5 ayant une longueur
L = 5 qui est soumise à un flux thermique caractérisé par un nombre de Rayleigh
de 1015. On retrouve sur cette figure trois courbes correspondant respectivement à la
solution estimée par l'équation (3.13) issue de l'approche analytique, la solution fournie
par Fluent et celle du modèle numérique à valider.

1 1 1 l S

/ /

0.8 -

/ /
/ /
S
ta o.6 - /
/
/
/
i 0.4
/ /

/ /
0.2 - / /
- » - modèle numérique
modèle Fluent
méthode intégrale
- » * ^ ' i
xlO4
0 1 2 3 4
Température à la paroi chaude Tw

FiG. 3.6 - Validation du modèle : T en fonction de y, cj> = 0.5, Ra = 1015, L = 5,

D = 0.001

On constate que la température obtenue grâce à la démarche analytique est quelque

peu surestimée, mais s'avère tout de même être une bonne indication de l'ordre de
 4]

10

Résultats Ra=1015
T = 3.62 10-4 y~° W,\
CD" Fluent Ra = 1015
Résultats Ra = 1016
cC
M Fluent Ra = 1016
o Résultats Ra = 10" ^"
■ r-l
7
S Fluent R a = 10'
^
cd ^
a T = 1.68 1 0 4 y 0 5 8
cd
i—*
-cd
d)
t*
T = 7.84 ÎO 5 y^-57
+J
Kl
?H

0)
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Position verticale dans le stratifié, y

FlG. 3.7 - Validation du modèle : T en fonction de y, 0.5,Ra= 1015-1017, L = 5

D = 0.001

grandeur des quantités en jeu. Les values obtenues grâce à l'équation 3.13 sont environ
10% supérieures aux résultats calculés d'après les approches numériques basées sur les
volumes finis. De plus, on remarque que les résultats de Fluent et du modèle numérique
coïncident fortement. D'ailleurs, cette similarité entre les deux modèles numériques est
particulièrement grande pour l'ensemble des nombres de Rayleigh traités, comme en
témoigne la figure 3.7 qui illustre la variation de température à la paroi chaude calculée
d'après les deux approches.

Les simulations indiquent bien que Tw évolue comme une puissance de y. Indépen­
damment de la porosité du milieu, pour le modèle numérique developé dans ce mémoire
l'exposant b moyen est 0.56 alors que dans Fluent ce paramètre correspond plutôt à
0.59, tel que montré aux tableaux 3.5 et 3.6.

Une grande importance est portée au calcul de la température maximale, car il s'agit
du critère à minimiser lors de l'optimisation. Par conséquent, il est primordial de bien
évaluer cette quantité lors de la modélisation. Le tableau 3.7 compare les valeurs de
températures maximales obtenues à la suite des deux types de simulation numérique.
Les différences relatives entre les résultats fournis par le calcul Fluent et ceux du modèle
 42

TAB. 3.5 - Coefficients de l'expression f(y) = ayb, (f> = 0.5, N = 1, L = 5,£> = 0.001
Fluent Modèle
Ra a b a b
14
10 7.9812E-4 0.617 7.9498E-4 0.615
15
10 3.6148E-4 0.589 3.6193E-4 0.586
10i6 1.6749E-4 0.587 1.6844E-4 0.583
1017 7.7679E-5 0.586 7.8418E-5 0.571
1018 3.6055E-5 0.586 3.6832E-4 0.548

TAB. 3.6 - Coefficients de l'expression fw(y) = ayb , (p = 0.90, N = 1, L = 5, D = 0.001

Fluent Modèle
Ra a b a b
1Q14
1.0942E-3 0.588 1.0913E-3 0.552
15
10 5.0787E-4 0.589 5.0939E-4 0.554
6
10i 2.3537E-4 0.587 2.3667E-4 0.552
1017 1.0916E-4 0.587 1.1004E-4 0.546
1018 5.0663E-5 0.588 5.1399E-5 0.521

de cette étude sont de l'ordre de 1% et moins. Cette concordance est aussi valable pour
de grandes porosités comme l'indique le tableau 3.8.

TAB. 3.7 - Validation du modèle numérique en considération de T max ,<p — 0.5, N = 1,

D = 0.001 . . .
Ra fmax Fluent fw Modèle % erreur
1014 8.1314E-4 8.0943E-4 0.45 %
15
10 3.6405E-4 3.6379E-4 0.07 %
6
10l 1.6879E-4 1.7009E-4 0.07 %
1017 7.8352E-5 7.8717E-5 0.04 %
1018 3.6367E-5 3.G899E-5 1.44 %

Le transfert de chaleur global extrait de la plaque chaude peut être quantifié par
la température moyenne à la paroi. Ce paramètre adimensionnel est calculé selon l'ex­
pression suivante :
 43

T A B . 3.8 - Validation du modèle numérique en considération de Tmax, 4> = 0.90, N = 1,

D = 0.001
Ra r m a x Fluent fw Modèle % erreur
14
10 1.1008E-3 1.101E-3 0.01 %
15
10 5.1093E-4 5.1328E-4 0.45 %
16
10 2.3719E-4 2.3846E-4 0.53 %
17
10 1.1010E-4 1.1086E-4 0.68 %
18
10 1.0890E-5 5.1814E-5 1.39 %

r m o î / = f1 fw dy (3.14)
Jo

où T moy représente la température moyenne à la paroi et Tw, la température à la paroi

chaude. Au tableau 3.9, on retrouve une liste des Tmoy obtenus pour différents Ra.
L'écart entre les deux simulations (Fluent et notre modèle), en regard du Tmoy, se situe
autour de 1 % et ce pour une vaste gamme de nombre de Rayleigh. Ce dernier indicateur
permet de confirmer la justesse des calculs et de valider le modèle numérique employé
dans ce mémoire.

T A B . 3.9 Validation du modèle numérique en considération de Tmoy, 0 = 0.90, N = 1,

1 = 5, D = 0.001
lia fmoy Fluent fmoy Modèle % erreur
10" 0.4458 0.4528 1.54 %
10ir> 0.9686 0.9778 0,94 %
16
LO 2.0943 2.0864 0.03 %
H)17 4.5202 4.4703 1.11 %
18
LO 9.7391 9.3288 4.39 %
 A4

TAB. 3.10 - Validation du modèle numérique en considération de Tmoy, 4> = 0.5, N = 1,

1 = 5, D = 0.001
Ra fmoy Fluent fmoy Modèle % erreur
10" 0.458 0.448 2,92 %
15
LO 0.978 0.968 0,94 %
1016 2.086 2.0943 0.38 %
17
10 l,4568E-3 1.4568E-3 7,56 %
1018 l,4568E-3 l,4568E-3 7,56 %

3.5 Conclusion

Le problème thermofluide que constitue l'échangeur poreux est résolu grâce à une
approche numérique basée sur les volumes finis. Au cours de ce chapitre, les principaux
éléments du modèle ont été présentés en détail. Dans un premier temps, on traite
du maillage et de la discrétistion du domaine. Puis, une analyse d'indépendance de
maillage est présentée. Les résultats du modèle numérique sont ensuite comparés à ceux
issus du logiciel commercial Fluent. La grande similitude, pour la plage d'opération
considérée, entre les résultats de Fluent et ceux de cette étude permet de valider le
modèle employé. De plus, les résultats numériques indiquent que les expressions issues
de la méthode intégrale, présentée à la section 2.3, permettent d'estimer la valeur de
la température maximale avec une précision d'environ 10%. Pour un empilement de
couches données, il est désormais possible de calculer numériquement avec certitude la
température maximale au sein du stratifié. L'étape suivante consiste à minimiser cette
température en optimisant l'architecture poreuse de l'échangeur. Cette optimisation est
réalisée grâce à un algorithme génétique.
 Chapitre 4

Algorithme génétique

Aux chapitres 2 et 3, les modèles mathématiques et numériques d'un échangeur de

chaleur formé d'un stratifié poreux ont été présentés. À ce stade-ci, les performances
thermiques du système peuvent donc être calculées en fonction de l'architecture du
refroidisseur. La prochaine étape consiste à optimiser la configuration des couches po­
reuses. Pour ce faire, l'approche préconisée est une technique de recherche s'appuyant
sur les algorithmes génétiques (AG). Au cours de ce chapitre, les caractéristiques de
fonctionnement de l'AG seront abordées soit : la codification d'une population (sec­
tion : 4.2.1), l'évaluation (section : 4.2.4) et le classement des individus (section : 4.2.5),
la reproduction (section : 4.2.11), la recherche locale (section : 4.2.8), la mutation (sec­
tion : 4.2.9) et finalement la convergence de la recherche (section : 4.2.11).

45
 46

4.1 Introduction

L'algorithme génétique (AG) est une technique de recherche qui s'inspire des prin­
cipes biologiques de l'évolution. L'AG agit sur une population de solutions potentielles
en appliquant des principes de sélection naturelle afin de générer un meilleur design.
Pour chacune des nouvelles générations créées, un ensemble d'individus est sélectionné
en fonction de leur niveau de performance par rapport au problème d'optimisation. Ces
individus, supérieurs aux autres, sont ensuite combinés en utilisant des opérateurs sem­
blables à ceux que l'on retrouve dans le monde de la génétique. Au fil des générations,
ce processus permet au groupe d'évoluer et de s'adapter à son environnement.

La figure 4.1 illustre les étapes définissant la procédure d'optimisation par algo­
rithme génétique. Tout d'abord, une population initiale est créée aléatoirement, chaque
individu représentant une structure particulière d'un stratifié, c'est-à-dire un design.
Des ajustements sont portés aux individus dont la configuration ne respecte pas le sens
physique du problème. Ensuite, ces individus sont évalués puis classés selon un indice de
performance. Parmi les meilleurs, certains sont sélectionnés pour former un ensemble
de parents. Grâce à un opérateur de reproduction, des enfants sont issus du groupe
de parents. Par la suite, des opérations de mutation et de recherche locale modifient
ces enfants de manière à bien couvrir l'ensemble du domaine de design. C'est ainsi
qu'une nouvelle population, possédant une partie des meilleures caractéristiques des
générations précédentes, est créée. Ce processus se répète jusqu'à ce qu'un critère de
convergence soit respecté.

La routine d'optimisation par algorithme génétique utilisée au cours de cet ouvrage

est construite à partir d'opérateurs issus d'un bloc de scripts libres développés dans
MATLAB [44].
 47

Population initiale

Décodage

Correction

Évaluation Nouvelle
^ population

Y
Sélection

t
Croisement

y
Recherche
locale

y
Mutation

y
Décodage

y
Correction
Meilleurs
individus y
Nouvelle
F population

>r
Non
( Convergence

"J~Oui

0
FiG. 4.1 - Schéma de fonctionnement de l'algorithme génétique
 48

4.2 Mécanismes de l'algorithme génétique

4.2.1 Génération de la population initiale

La première étape consiste à construire une population initiale constituée de Nina

individus. Chaque individu représente une architecture particulière d'un stratifié, et
à l'image du monde de la génétique, ces caractéristiques sont codées sous la forme
de chromosomes. Une constitution spécifique, nommée phenotype, correspond à un
arrangement de gènes, ou degrés de liberté. Chaque individu possède un nombre de
gènes égale à deux fois le nombre de couches TV du stratifié, un gène correspond au
matériau et un second correspond à la porosité de la couche. Le regroupement de gènes
reliés aux matériaux se nomme chromosome de matériaux et similairement, l'ensemble
des gènes de porosité se nomme chromosome de porosité. Ces gènes sont encodés grâce
à une suite de nombres binaires. Par conséquent, l'alphabet des gènes comprend deux
valeurs, soit 0 soit 1. Le domaine associé à une variable est donc fractionné en 2Nbit,
où Nbn représente le nombre de bits associé à un gène. Dans ce mémoire, la porosité
et les matériaux sont considérés comme des variables discrètes codées sur 5 bits et 2
bits respectivement. La figure 4.2 illustre la composition d'un phenotype désignant un
stratifé ayant deux couches (N = 2). Les premiers individus sont générés seulement
d'après la longueur totale des chromosomes Lin(i sans égard au sens physique que peut
représenter la suite binaire. Les gènes de départ sont donc sélectionnés aléatoirement.
Par exemple, l'individu illustré à la figure 4.2 possède une valeur de Linc[ égale à 14.
Fait à noter, la longueur physique d'un stratifié n'est pas directement traitée sous
la forme d'une variable. Cette dernière est plutôt préalablement fixée à une certaine
valeur. Lors de l'optimisation, l'algorithme génétique est libre de dimininuer la taille du
stratifié en imposant des couches vides. Ce mécanisme permet d'optimiser indirectement
la longueur du stratifié soumis à certaines contraintes. Une fois que l'ensemble des
variables est défini, l'étape suivante consiste à décoder les chromosomes afin de connaître
l'architecture poreuse des différents individus.
 49

individu
^
chromosome c h r o m o s o m e de
de m a t é r i a u porosité
~>»/
gène de gène de
matériau porosité
bit

0 1 0 0 0 0 0 10 0 1 1 1 0
couche1 couche couche: couche,

FlG. 4.2 - Représentation binaire d'un individu, N = 2

4.2.2 Décodage

Sans indication sur la structure et la signification des chaînes de bits, on ne peut

pas associer une caractéristique physique du stratifié telle une porosité à un gène. Par
conséquent, pour décoder un chromosome on doit connaître le domaine et le système
numéral de chaque variable. Cette étape s'effectue à la suite de la création d'une nouvelle
population, comme l'indique la figure 4.1.

Comme il a été mentionné à la section 4.2.1, la porosité, qui peut prendre une valeur
comprise entre 0 et 1, est codée sur 5 bits et est traitée comme une variable discrète.
Le fractionnement du domaine est linéaire ce qui implique que la différence entre deux
valeurs de porosité successives est constante. Ainsi, la précision sur cette variable est
fonction du nombre de bits utilisé pour coder un gène. Dans ce cas, la résolution est de
1/25, soit environ 3%. Le système numéral employé pour représenter ce type de gène est
le code Gray [45]. Ce choix est particulièrement approprié aux algorithmes génétiques,
car les nombres d'un même voisinage diffèrent très peu dans leur représentation binaire.
Cette représentation permet d'améliorer la convergence et le temps de calcul de l'AG
soumis à une optimisation globale sur des paramètres continus [46].

Les gènes reliés aux matériaux sont encodés de la même manière que les gènes de
porosité. Ainsi, grâce au système numéral employé, deux variables correspondant à des
matériaux positionnés successivement dans la liste 2.1 ne diffèrent que par un seul bit.
 50

Étant donné que les gènes de matériaux sont codés sur deux bits, quatre combinaisons
sont possibles. La quatrième combinaison est réservée à une couche vide (cfi = 1).
Les matériaux un à trois correspondent respectivement à de l'aluminium, du cuivre
et de l'acier. La figure suivante, qui reprend l'exemple présenté à la section précédente
(figure 4.2), illustre la manière dont on décode l'architecture d'un stratifié à partir d'une
chaîne binaire.

matériau porosité
~>s~~
couche. couche, couche couche,

0 1 0 0 0 0 0 10 0 1 1 1 0

Représentation en code Gray

matériaux porosité
code Gray correspondance code Gray correspondance

0 0 -> 1 -+• aluminum 0 0 0 0 0 - ► 1 -*• 1/32 = 0.0313

0 1 -> 2 - ► cuivre 0 0 0 0 1 —> 2 - ► 2/32 = 0.0625
11 -► 3 -► fer 0 0 0 1 1 -> 3 -+> 3/32 = 0.0938
10 - 4 - couche vide 0 0 0 1 0 - * 4 -*" 4/32 = 0.1250

1 0001 _► 31 _► 31/32 = 0.9688

1 0000 32 32/32 = 1.0000

Représentation physique

eu al 0.1250 0.3750
t t t t
matériau ï matériau 2 porosité, porosité 2

FlG. 4.3 - Exemple de décodage d'une chaîne binaire, N

Selon cette modélisation, un stratifié constitué de dix couches (N=10) peut acquérir
70
2 architectures différentes. Dans ce contexte, l'emploi d'un AG s'avère être un outil
puissant pour trouver efficacement une solution à un problème d'optimisation.
 r>i
4.2.3 Correction des individus

Une fois le décodage des chaînes binaires terminé, il est possible que l'architecture
d'un ou plusieurs individus ne respecte pas le modèle physique. Par exemple, une couche
vide doit être positionnée à l'extrémité de la structure pour que la résolution numérique
soit valide. Si tel est le cas, des corrections doivent être apportées à l'architecture
poreuse afin que le calcul thermofluide soit significatif. Un nouvel opérateur a donc
été créé afin de remplir cette tâche. Ce dernier inspecte chaque nouveau phénotype
afin de déceler de potentielles anomalies. Deux types d'ajustement sont apportés aux
individus. Premièrement, le dernier matériau de la liste désigne une couche vide dont
la porosité correspondante est <>
/ = 1. Par le fruit du hasard, lors de la génération de
la population initiale, ou suite à une mutation (voir section : 4.2.9), il se peut que la
quatrième combinaison de matériau soit attribuée à une couche ayant une porosité autre
que 100%. Dans ce cas, un nouveau matériau, sélectionné aléatoirement, est assigné à
la couche et la porosité est conservée. Dans un même ordre d'idées, lorsqu'une variable
de porosité atteint la valeur maximale (<j) = 1), le gène de matériau couplé avec cette
dernière doit correspondre à une couche vide. Par conséquent, une strate vide apparaît
uniquement lorsqu'une des variables de porosité l'impose.

Deuxièmement, un stratifié doté d'une couche vide au cœur de sa structure est dé­
pourvu de sens physique. En effet, une strate dont la porosité est de 100% fait office
de frontière adiabatique, car la conductivité thermique de cette dernière est nulle. Par
conséquent, les couches qui succèdent une couche vide ne contribuent plus au transfert
de chaleur et deviennent alors thermiquement inutiles. Lors de la résolution thermo­
fluide, ce type de configuration ne pose pas de problème. Cependant, l'évaluation de
la masse et du coût est biaisée par la contribution de la portion superflue du stratifié.
Également, l'AG associe une valeur de performance à un phénotype, sans connaître
la nature du problème traité (voir section 4.2.4). Lors d'un problème non contraint,
l'algorithme génétique tentera alors d'optimiser cette portion en vain, ce qui aura pour
effet d'accroître le nombre d'évaluations d'individus inutilement. C'est pourquoi un
opérateur spécifique à ce problème a été conçu. La stratégie est simple : une fois que
les ajustements relatifs au couplage des variables sont réalisés, l'opérateur parcourt les
chromosomes à la recherche de couches vides. Lorsqu'une couche ayant une variable
4> = 1 est rencontrée, l'opérateur élimine les couches subséquentes en remplaçant les
gènes de matériaux et de porosité par ceux associés au vide. En d'autres termes, le
stratifié est tronqué à partir de la strate vide la plus près de la paroi chaude. Cette
 52

correction permet d'éliminer rapidement les couches considérées superflues. Ce travail

aurait pu être réalisé par le biais de la mutation et du croisement, mais l'ajustement
s'étalerait possiblement sur plusieurs générations ce qui induirait un coût de calcul
important. Indirectement, grâce à ce mécanisme de correction, l'algorithme génétique
optimise la longueur du stratifié. Cependant, le processus d'optimisation ne permet pas
d'allonger l'échangeur au-delà de la valeur maximale fixée d'après les conditions initiales
du problème. La figure 4.4 présente un exemple de corrections apportées aux individus.
Le phénotype initial, que l'on retrouve au haut de la figure 4.4, comporte certaines
anomalies qui doivent être corrigées avant que l'individu ne soit évalué. Premièrement,
la configuration poreuse comprend une porosité de 100% à la troisième couche alors
que le matériau correspondant ne désigne pas une couche vide (matériau 4). Le gène
de matériau de la troisième couche est donc replacé pour respecter la règle de couplage
des variables. De plus, d'après le chromosome de matériau, la quatrième couche devrait
être vide. Cependant, la porosité de cette dernière n'est pas égale à 1. Alors, le gène
de matériau de la quatrième couche est remplacé pour désigner aléatoirement un des
trois métaux conducteurs de la liste 2.1. Cette dernière étape est illustrée à la deuxième
ligne de la figure 4.4. Finalement, le stratifié est tronqué à partir de la première couche
vide rencontrée en s'éloignant de la paroi chaude. Les strates 3 et 4 sont par conséquent
substituées par des couches vides.

Chromosome Chromosome de
de matériau porosité

1 3 2 4 3 0.1 0.7 1 0.0 0.5

Phénotype initial

1 3 4 1 3 0.1 0.7 1 0.0 0.5

Couplage des variables

1 3 4 4 4 0.1 0.7 1 1 1

Tronquage de la séquence

FlG. 4.4 - Exemple de corrections apportées à un individu, N = 5

4.2.4 Evaluation des individus

Après que les corrections aient été apportées à la population, les individus sont prêts
à être évalués. Le processus d'optimisation repose principalement sur la performance
des individus, désignée par une fonction objectif symbolisée par f(x) où x représente
un individu. L'algorithme génétique fonctionne donc sans connaître la nature du pro­
blème. En effet, l'AG traite seulement un ensemble de solutions auxquelles des valeurs
numériques de f(x) sont associées. La qualité d'une solution est donc mesurée par la
valeur de f{x). L'algorithme minimise la valeur de la fonction objectif. Ainsi plus f(x)
est faible, meilleur est considéré le phénotype.

Chaque phénotype correspond à une solution différente au problème thermofluide.

La performance des individus se calcule à partir de propriétés issues de la résolution du
modèle numérique présenté au chapitre 3. L'algorithme génétique a pour but de déceler
l'architecture poreuse optimale qui minimise la température maximale à la paroi. La
propriété qui fait gage d'un transfert thermique efficace est la température maximale
dans l'échangeur. Par conséquent, lors d'une optimisation sans contrainte, les stratifiés
caractérisés par une faible valeur de Tmax s'avèrent être les plus performants. La qualité
de ces solutions potentielles en rapport à leur environnement se quantifie alors par la
fonction objectif f(x) suivante :

f(x) = fmax (4.1)

La formulation précédente de la fonction objectif (4.1) ne tient pas compte de la

masse et du coût du stratifié dans le calcul de la performance. Dans le cas de problèmes
contraints, on doit inclure les contraintes à même la fonction objectif. Pour ce faire,
une technique de pénalisation est employée. La fonction de pénalisation altère la per­
formance des individus qui violent la ou les contraintes en augmentant la valeur de la
fonction objectif. Ceci fait en sorte que ces derniers sont défavorisés. La pénalité aura
alors pour effet de limiter les chances de reproduction d'un individu fautif, de telle sorte
que la configuration de ce dernier tendra éventuellement à disparaître de la population
au fil des générations (voir section 4.2.7). Lorsque la masse et/ou le coût associés à
une structure excèdent certaines limites prescrites, on augmente la valeur de la fonction
objectif et l'équation 4.1 devient :
 54

f(x) = fmax \l + em m a x ( ( - = - l),0) f ec max ( l),0) (4.2)

,MC0

où eTO et ec sont des paramètres qui accentuent l'effet de la pénalité. Mco et Cco re­
présentent respectivement les contraintes de masse et de coût. Lorsque les contraintes
ne sont pas respectées ( -^- > 1 et/ou ■£- > 1 ) la valeur de T m a x est multipliée
par un scalaire proportionnel au carré de l'écart entre les propriétés du stratifié et les
propriétés maximales admises. Cette opération a pour but de réduire les probabilités
que ce type de configuration se reproduise au sein de la population. Dans le cas où M
et C demeurent sous les limites prescrites, la pénalité disparaît et la fonction objectif
correspond uniquement à T m a x .

Les paramètres em et ec sont ajustés de manière à éliminer efficacement les solutions

qui excèdent fortement les limites imposées. Cependant, l'algorithme génétique doit
tout de même être en mesure d'inspecter le domaine situé près des contraintes, car
il est possible qu'une solution optimale émerge de l'amélioration d'une configuration
faiblement pénalisée. Des valeurs trop élevées de em ou ec ont pour effet de trop limiter
la recherche dans l'espace de design.

4.2.5 Classement des individus

Après avoir évalué les individus, on classe ces derniers selon leur qualité. Cette
étape a pour but de classer les individus selon leur performance par rapport aux autres
membres du groupe. Ce classement sera ensuite employé lors de la sélection qui est dé­
crite à la section suivante. Généralement, la valeur numérique calculée à l'équation (4.2)
est utilisée seulement dans une étape intermédiaire servant à déterminer la performance
absolue des individus d'une population par rapport à l'environnement. Une autre fonc­
tion, nommée fonction de classement, sert à transformer les valeurs de la fonction
objectif f(x) en mesure de classement relatif :

F(x)=g(f(x)) (4.3)
 55

où / (x) désigne la fonction objectif , g (x), la fonction de classement et F (x) représente

la performance relative. Globalement la valeur de la fonction de classement indique
le nombre d'enfants qu'un individu peut espérer produire au cours de la prochaine
génération. Dans la littérature, on retrouve plusieurs types de transformations. Une
manière simple d'exprimer la fonction de classement est :

Fi-ù = ^S:\ ; (4-4)

où Nind correspond à la taille de la population et Xi représente le phénotype d'un

individu. Les transformations linéaires comme celle énoncée à l'équation (4.4) sont sus­
ceptibles de provoquer une convergence locale rapide, mais souvent prématurée [47].
En effet, le processus de sélection, décrit à la section 4.2.7, repose sur la valeur de
performance. Dans ces circonstances, le nombre d'enfants issus d'un individu est pro­
portionnel à la valeur de F(x). S'il n'y a pas de limitation sur la plage de valeur de
performance, un phénotype très performant qui apparaît au cours des premières généra­
tions peut dominer la reproduction et causer une convergence hâtive vers un minimum
local. De plus, selon cette approche, il s'avère difficile de déceler le meilleur individu
d'une population dont les membres possèdent des valeurs de performance similaires.
Afin de régler ce problème potentiel, au cours de cette étude, une variante de l'équa­
tion (4.4) est utilisée (éq. (4.5)). Le classement, au lieu de simplement se calculer grâce
à la valeur de performance relative, se détermine plutôt d'après le rang d'un individu,
après que la population ait été triée en ordre décroissant de valeur de fonction objectif
/ (x) [48]. La transformation employée se définit comme :

fw-4^ («)
La qualité des différentes solutions proposée est désormais quantifiée par le résultat
de l'expression 4.5. Une fois le classement achevé, l'étape suivante consiste à sélectionner
les phénotypes aux fins de la reproduction à partir de la liste qui vient d'être construite.
 56

4.2.6 Sélection

En biologie, les organismes évoluent au fil des générations afin de mieux s'adap­
ter aux changements de leur environnement. Les mécanismes à la base de l'évolution
sont la sélection naturelle et la variation génétique. Face à des changements de leur
environnement, certains individus d'une population possèdent des attributs génétiques
qui favorisent leur survie, et du même coup, leur chance de se reproduire. Ces derniers
transmettent alors leur bagage génétique à leurs enfants qui forment la nouvelle généra­
tion. En répétant ce processus de sélection sur plusieurs générations, la population tend
à s'améliorer. C'est sur ce principe fondamental de la nature que se base la recherche
de l'algorithme génétique.

Avec l'AG, la sélection consiste à déterminer le nombre de fois qu'un individu est
choisi pour se reproduire, et par le fait même, le nombre d'enfants que cet individu
produira. Cette opération est complémentaire au classement des individus impliqués
à la section précédente. Le mécanisme de la sélection se base sur l'ordonnancement
des solutions formulées à l'équation 4.5 pour établir les probabilités de reproduction
des membres de la population. Ainsi, l'individu qui détient la plus grande valeur de
F(x) possède les meilleures chances de se reproduire. Toutefois, selon la technique de
sélection employée, il n'est pas automatiquement assuré que ce dernier sera choisi.

Plusieurs algorithmes de sélection ont été développés dans la littérature. Parmi les
plus performants, on retrouve la sélection stochastique universelle [47] (SUS). Cette
dernière est employée au cours de cette étude. La méthode SUS utilise les valeurs de
performance relative F(x) issues du classement, calculées grâce à l'équation (4.5), pour
déterminer la probabilité de reproduction des individus. Chaque individu est positionné
sur un segment de longueur S, où S correspond à la somme des performances relatives,
soit S = Yld=\dFi (x). L'espace qu'occupe un individu Sj sur ce segment est propor­
tionnel El SB. performance au sein du groupe tel que :

F
s = m (46)

Afin de mieux visualiser ce processus, la figure 4.5 schématise l'utilisation de l'al­

gorithme de sélection stochastique universelle. On dispose les phénotypes en ordre dé-
 57

croissant de longueur. Finalement, pour désigner quels individus seront choisis, cette
méthode utilise Nsei pointeurs simultanément, où Nsei correspond au nombre de sélec­
tion requise. Ces pointeurs sont également espacés entre-eux d'une distance équivalente
à S/Nsei. Le premier pointeur se positionne aléatoirement dans l'intervalle [0,S/Nsei\.
Cette stratégie permet de réduire le risque qu'un même phénotype soit sélectionné plu­
sieurs fois comme c'est le cas avec d'autres méthodes de sélection comme la roulette.
Grâce à la répartition des probabilités selon la performance relative, les individus mieux
adaptés à leur environnement sont naturellement avantagés et donc plus susceptibles
de contribuer à la reproduction. Toutefois, il est possible que certains phénotypes de
moindre qualité puissent être sélectionnés lors du positionnement des pointeurs.

position aléatoire

Y s, =0.33 s 2 =0.21 s3 = 0.16 s, = 0.12 s 5 =0.11 s6= 0.07

1 2 3 4 5 I 6
0.0 0.33 0.54 0.70 0.82 0.93 1.0

X ™ x *± X - î
pointeur 1 pointeur 2 pointeur 3 pointeur 4

FlG. 4.5 - Sélection des individus par la méthode SUS, Nind=6 et iVse/=4

4.2.7 Reproduction

L'opérateur à la base de la production de nouveaux chromosomes dans l'algorithme

génétique est le croisement. À l'instar du monde vivant, la reproduction génère de
nouveaux individus qui possèdent une partie du bagage génétique des deux parents.
Comme pour la plupart des éléments de l'AG, la mécanique du croisement contient
une grande part d'aléatoire. De plus, il existe une vaste gamme d'opérateur de croise­
ment. En effet, on retrouve des opérateurs spécifiques à certains types de gènes ainsi
que pour différentes applications. Pour des variables discrètes codées sous une forme
binaire, on choisit un croisement multipoint. Cette technique de reproduction, illustrée
à la figure 4.6 est adoptée dans cette étude. Tout d'abord, on sélectionne aléatoire­
ment un certain nombre de positions de croisement sans redondance dans l'ensemble
[1,2,..., Lind — 1], où Lin(i représente le nombre de bits que contient un individu (voir
section 4.2.1). Le nombre de points de croisement est désigné par la variable m. Ces
points sont ensuite triés en ordre ascendant. Les bits contenus entre deux points de
 58

croisement successifs sont échangés entre les parents pour produire les enfants. La sec­
tion entre le début de la chaîne binaire et le premier croisement demeure inchangée.
Ce processus est illustré à la figure 4.6. On peut constater que la coupure de la chaîne
binaire se fait au gré du hasard sans tenir compte de la longueur des gènes. Ainsi, la
nature discontinue de cette stratégie permet d'encourager l'exploration du domaine de
design plutôt que de favoriser une convergence vers un individu hautement performant
au début d'une recherche.

eu al <() = 0.125 <|> = 0.375

parent 1 0 1 0 l 0 0 0 0 1| 0 0 1 1 1

lie eu 0 = 0.813 $ = 0.531

parent 2 1 1 0 ,1 1 0 1 1,1 1 1 0 0

1
eu eu 0 = 0.844 0 = 0.313
enfant 1 0 1 0 1 10 1 1 0 0 1 1 1 1

fe al 0 = 0.625 <)> = 0.500

enfant 2 1 1 0 0 0 0 0 11 1 1 0 0 0

FlG. 4.6 - Croisement de deux individus par la méthode mufti-points, m = 3

4.2.8 Recherche locale

Dans le but d'améliorer la recherche de l'optimum par l'algorithme génétique, on

implante une routine d'amélioration locale. Les opérateurs nécessaires pour effectuer
cette étape ne figurent pas parmi le groupe de scripts libres développés dans MAT-
LAB [44]. De nouvelles fonctions ont donc été créées spécialement pour réaliser cette
tâche. L'hybridation de l'AG permet de déceler des minimums locaux autour des in­
dividus d'une population, ce qui bonifie globalement leurs qualités. L'AG fonctionne
sur des principes probabilistes. Avec l'AG, il est donc parfois fastidieux de détecter,
uniquement à l'aide des opérateurs de croisement et de mutation, un optimum qui se
situe dans l'entourage immédiat d'un individu. L'hybridation a pour effet de pallier à
ce problème et de réduire le temps de calcul [32].

La recherche locale s'applique sur l'ensemble des individus d'une population, et ce,
périodiquement à toutes les Niocaie générations. La recherche locale ne s'effectue pas
à chaque génération en raison du temps de calcul associé à cette opération. De plus,
cette étape s'avère particulièrement utile après plusieurs générations, à l'approche d'un
optimum. La stratégie consiste à évaluer les voisins des individus dans l'espace de
design afin de vérifier la possibilité d'améliorer la fonction objectif f(x). Le voisinage
se définit à partir du phénotype d'un individu auquel on applique de légères variations.
Si ces modifications conduisent à une diminution de f(x), la structure initiale est alors
remplacée par le résultat de la recherche locale.

La représentation génétique de l'échangeur de chaleur comprend deux chromosomes ;

l'un représente les matériaux de l'empilement, l'autre correspond à la porosité des dif­
férentes couches. Le voisinage de ces deux types de variables est traité différemment.
Dans un premier temps, on génère une série de voisins à partir du chromosome de
porosité en gardant celui des matériaux inchangés. Pour chaque individu, on tire aléa­
toirement un chiffre x compris entre 1 et TV (nombre de couches du stratifié). Ce chiffre
sert à déterminer combien de gènes seront inspectés. Pour chacun des gènes sélection­
nés au hasard, on génère deux voisins à partir des valeurs supérieures et inférieures
de la variable associée à la caractéristique en question. Dans le cas où le gène choisi
correspond à une frontière du domaine, on sélectionne les deux valeurs qui précèdent
ou qui succèdent à la variable. À la figure 4.7, un schéma du domaine de la recherche
locale est présenté. Quatre voisins sont créés à partir de l'inspection des deux variables
de porosité associée à deux couches. Selon cette approche, si tous les gènes de porosité
sont choisis, un maximum de 2N voisins peuvent être générés. Par la suite, ces nouveaux
phénotypes sont évalués. Le meilleur des nouveaux designs est comparé au phénotype
initial et si la modification est avantageuse alors la structure génétique de l'individu
est remplacée par celle de son voisin. Ce procédé se répète pour chacun des membres
de la population. Un exemple de recherche locale appliquée à la porosité est illustré à
la figure 4.8. Etant donné que l'on détermine aléatoirement la taille du voisinage, le
temps de calcul associé à l'hybridation de l'AG n'est pas constant d'une génération à
une autre.
 60

4>2+A(|>2
V
4>, - A<t>, î <!), + A * ,

I
o At
<t>2 - >2

<K
FiG. 4.7 - Schéma du domaine de la recherche locale appliquée au chromosome de
porosités, N = 2, x = 2

Chromosome Chromosome
de matériau de porosité
-^ ^
^

1 3 2 1 3 0.1 0.4 0.8 0.1 0.5

Phenotype initial

1 3 2 1 3
t 0.1 0.4 0.9 0.1 0.5
Voisin supérieur

1 3 2 4 4 0.1 0.4 0.7 0.1 0.5

Voisin inférieur

FiG. 4.8 - Exemple de recherche locale appliquée au chromosome de porosité, N = 6,

x= 1
 ()l

Dans un deuxième temps, on construit le voisinage associé au chromosome des ma­

tériaux sans toucher à celui de la porosité. L'approche employée est différente de la
précédente. Pour chaque individu, la stratégie consiste d'abord à regrouper les strates
faites du même matériau. Ensuite, les voisins sont générés à partir des différentes com­
binaisons possibles d'agencement de matériaux auxquelles on greffe le chromosome de
porosité. Ainsi, les nouveaux phénotypes créés possèdent tous la même structure po­
reuse, car le chromosome d'origine demeure intact. Selon cette approche, un maximum
de six voisins par individus peuvent être créés, en condésirant le tableau 2.1 qui contient
trois types de matériau. La meilleure de ces solutions potentielles est comparée au phé-
notype de base et s'il y a amélioration, l'individu est remplacé par son voisin plus
performant. Un exemple de recherche locale effectuée sur le chromosome des matériaux
est présenté à la figure 4.9.

4.2.9 Mutation

Dans la nature, la mutation est considérée comme une des forces dominantes de
l'évolution. Ce processus modifie aléatoirement le matériel génétique d'un organisme, ce
qui a pour effet de créer une nouvelle structure génétique. Certaines mutations peuvent
engendrer des modifications qui altèrent la performance d'un individu. Par exemple,
considérons la mutation d'un gène de porosité qui aurait pour effet de réduire fortement
la porosité d'une des premières couches de l'échangeur. Cela aurait pour conséquence
de limiter le passage du fluide, provoquant une concentration de la chaleur à la paroi
chaude. Ainsi, la modification impliquerait une diminution de performance pour l'in­
dividu. Grâce à la sélection naturelle, les mutations défavorables sont graduellement
éliminées. En revanche, ces substitutions aléatoires peuvent parfois être fort bénéfiques.
Si tel est le cas, les modifications tendent à se propager au sein de la population au
fil des générations. Prenons le cas où une mutation change une couche d'acier en une
couche de cuivre. Le nouveau matériau étant bien plus conducteur que le précédent,
la répartition de chaleur au sein de l'échangeur sera améliorée. De ce fait, la perfor­
mance du design sera bonifiée, et le cuivre tendra à s'imposer au cours des générations
subséquentes.

Dans l'algorithme génétique, la mutation est nécessaire afin de bien diversifier la

recherche. Ce mécanisme peut produire des transformations de la structure qui ne
 (>2

Chromosome Chromosome
de matériau de porosité
-^ ^
^
r^
1 2 1 0.1 0.4 0.8 0.1 0.5

Phenotype initial

1 1 2 13 a, 0.4 0.8 0.1 0.5

Première combinaison

1 1 3 3 2 0.1 0.4 0.8 0.1 0.5

Deuxième combinaison

2 1 1 0.1 0.4 0.8 0.1 0.5

Troisième combinaison

- 13 1 3 1 1 1 0.1 0.4 0.8 0.1 0.5

Quatrième combinaison

Q
1 2 0.1 0.4 0.8 0.1 0.5

■
1

Cinqi îième combinaison

3 2 1 1 0.1 0.4 0.8 0.1 0.5

*
Sixième combinaison

FiG. 4.9 - Exemple de recherche locale appliquée au chromosome de matériaux, N = 5

 63

peuvent pas être obtenues par croisement. La mutation garantit que la probabilité de
générer une séquence spécifique n'est jamais nulle.

Plus concrètement, une fois que les enfants ont été générés, un opérateur parcourt
la chaîne binaire associciée à un individu et permute des bits aléatoirement selon une
certaine probabilité. Toutefois, les chances qu'un bit subisse une modification demeurent
faibles soit entre 0.001 et 0.01. Un taux de mutation trop élevé aurait pour effet de
faire diverger la recherche en raison de l'importance de l'altération des structures. La
figure 4.10 présente un exemple de l'impact d'une mutation sur un phénotype.

CU al (|> = 0.125 <j) = 0.375

0 1 0 0 0 0 0 10 0 1 1 1 0
avant mutation
Y
CU CU 0 = 0.125 $ = 0.625
0 1 0 1 0 0 0 10 1 1 1 1 0
après mutation

FlG. 4.10 - Exemple de mutation, N — 2

Fait à noter, la mutation peut potentiellement modifier le phénotype et conduire

à un design d'échangeur qui ne respecte pas le modèle physique. Par conséquent, à la
suite de cette étape, on doit ensuite procéder à une correction telle que présentée à la
section 4.2.3.

4.2.10 Nouvelle génération

Après avoir été soumis à la mutation, les enfants sont alors prêts à être réinsérés
dans la prochaine génération. Au cours de cette étude, la technique employée pour
construire la nouvelle population se fonde sur une stratégie multiélitiste. En effet, un
groupe constitué des Nm meilleurs individus issus de la génération précédente forme la
base de la génération en cours de construction. La qualité des individus sélectionnés
se détermine grâce à la fonction de performance relative F (x) (éq. 4.5). Ensuite, on
complète la population, contenant Ninci individus, par les enfants issus du croisement.
Le nombre d'individus constituant la progéniture est alors égale k Ne — Ninct — Nm,
où Ne désigne le nombre d'enfants. Lors de l'étape de la reproduction, Ne nouveaux
 M

phénotypes sont produits et la totalité de ces individus sont implantés dans la nouvelle
génération, et ce, sans tenir compte de la performance de ces derniers. De plus, la
proportion d'enfants par rapport aux élites doit être élevée afin d'explorer plusieurs
nouveaux designs et ainsi accentuer la diversité. Le rapport Ne/Nin^ utilisé dans ces
travaux est de 95%. Conséquemment, le nombre d'individus évalués nbevai au cours
d'une optimisation, sans l'application de recherche locale (voir section 4.2.8), est :

nbeval = Nind + (G - 1) Ne (4.7)

où G représente le nombre de générations totales effectuées. Pour un problème relati­

vement complexe le nombre de générations requises pour atteindre la convergence peut
parfois être grand (voir section 4.2.11). En outre, le temps nécessaire pour l'évaluation
d'un seul individu peut aussi s'avérer être élevé, notamment lorsque le nombre de Ray-
leigh est important (voir section 3.3). Dans le but d'optimiser le temps de calcul, une
base de données contenant l'historique de l'ensemble des phénotypes traités au cours
d'une même recherche est construite. Ainsi, si par le fruit du hasard un design déjà
évalué se retrouve dans la nouvelle génération, ses propriétés sont aussitôt récupérées
de la base données, ce qui évite d'effectuer le calcul thermofluide une deuxième fois.
Cette stratégie s'avère particulièrement efficace lorsque l'AG converge vers un minimum
global et qu'alors plusieurs individus partagent un matériel génétique identique.

4.2.11 Critère d'arrêt et convergence

La recherche se termine lorsque le critère de convergence est atteint. Ce dernier peut

prendre différentes formes. On peut tout simplement fixer un nombre de générations
maximales après quoi on interrompt le calcul. Si la valeur optimale de la fonction objectif
f(x) est connue, il est également possible d'arrêter l'optimisation lorsque l'on obtient
un résultat suffisamment près de l'optimum désiré. Une autre approche, fréquemment
employée, consiste à fixer le critère de converge à un nombre minimal de générations
consécutives sans amélioration du meilleur individu. En d'autres termes, l'algorithme
génétique poursuit sa recherche tant que la température maximale demeure inchangée
pendant plusieurs générations successives, ce qui constitue un critère de convergence.
Cette dernière approche est employée au cours de cette étude. Le nombre minimal de
 65

générations sans changement choisi est 250. Toutefois, un nombre total de générations
est aussi utilisé en guise de critère d'arrêt de sûreté.

4.3 Conclusion

L'optimisation d'un stratifié poreux employé comme refroidisseur en convection na­

turelle se révèle être un problème fort complexe. Le nombre de configurations possibles
s'accroît rapidement en fonction du nombre de couches. L'emploi d'une méthode de
recherche par algorithme évolutif s'avère donc nécessaire afin de déceler l'architecture
optimale dans un délai respectable. Au cours de ce chapitre, les principaux éléments du
fonctionnement de l'AG ont été abordés. Des opérateurs conventionnels comme le déco­
dage, la sélection, croisement et la mutation y ont été présentés. De plus, des fonctions
créées spécifiquement pour le problème thermofluide de ce mémoire ont été décrites
comme la correction, l'évaluation et la recherche locale. Une fois que l'algorithme géné­
tique est bien implanté, il est alors possible d'optimiser l'échangeur sous des contraintes
de masse et de coût. Les résultats de l'optimisation feront l'objet du prochain chapitre.
 Chapitre 5

Résultats de l'optimisation

Ce chapitre est consacré aux résultats obtenus suite à l'optimisation d'un stratifié
poreux à l'aide de l'algorithme génétique présenté au chapitre 4. La recherche d'un stra­
tifié qui minimise la température maximale à la paroi est d'abord réalisée sans imposer
de contraintes. Lors du processus, l'AG est donc libre de choisir parmi l'ensemble des
matériaux proposés au tableau 2.1 et la gamme de porosités, sans considération à la
masse et au coût du stratifié. Ensuite, l'optimisation est effectuée avec une contrainte
de masse, puis avec une contrainte sur le coût. Finalement, la température critique du
système est minimisée en considérant une combinaison de contraintes sur la masse et
sur le coût de l'échangeur.

66
 67

5.1 Optimisation sans contraintes

Tout d'abord, l'optimisation de la structure poreuse par rapport à la température

maximale T m a x est réalisée sans contraintes. Les paramètres em et ec de l'équation (4.2)
sont alors nuls. Par conséquent, la qualité d'un individu est uniquement fonction de la
répartition de température dans l'échangeur. Ainsi, une grande performance est attri­
buée à une structure ayant une faible température maximale, comme l'indique l'équa­
tion (4.1). Les variables du système sont la porosité et le matériau constituant les
couches poreuses. Le nombre de Rayleigh est fixé de même que la longueur du stratifié.
Les simulations sont réalisées avec un nombre de couche N suffisament grand de sorte
que l'architecture optimale émergeant de l'optimisation soit finement détaillée. De ce
fait, chaque stratifié possède dix couches d'égale épaisseur 1 .

Etant donné que l'optimisation par algorithme génétique est un processus stochas­
tique, le nombre de générations requises pour atteindre la convergence peut varier d'un
calcul à un autre, tout dépendant du parcours qu'effectue l'AG dans le domaine de de­
sign lors de la recherche. Ainsi, afin de s'assurer de la reproductibilité des résultats, cinq
recherches en utilisant les mêmes paramètres sont réalisées pour chaque cas considéré.
Par exemple, le design optimal 2 présenté au tableau 5.1 a été obtenu à cinq reprises
lors d'optimisations différentes nécessitant respectivement 707, 392, 479, 442 et 425
générations.

5.1.1 Influence de la longueur sur la résistance thermique

Cette section se consacre à l'étude de l'influence de la taille du stratifié sur les

performances thermiques. La dimension de l'échangeur varie grâce au paramètre L qui
désigne la longueur totale des couches. Lors des simulations, le nombre de couche du
stratifié est fixé à N = 10. De plus, l'ensemble des matériaux proposés au tableau 2.1
sont disponibles. Cependant, l'optimisation sans contraintes converge invénitablement
vers un échangeur constitué entièrement de cuivre, car, bien que ce matériau soit le plus
cher et le plus massif parmi les constituants possibles (tab. 2.1), il détient la meilleure
conductivité thermique.
'Un exemple d'application dimensionnel est présenté en annexe A
2
La première couche présentée dans le tableau 5.1 est adjacente à la source de chaleur
 68

Afin de mesurer l'effet de la longueur du stratifié sur la température maximale à la

paroi chaude T m a x , plusieurs optimisations sont réalisées pour différentes valeurs de L,
et ce, pour des nombres de Rayleigh variant de 1014 à 10 18 . Les résultats de cet exercice
sont résumés à la figure 5.1 où est illustrée la variation de la température maximale en
fonction de la longueur du stratifié. Après l'analyse de ces résultats, on constate que la
longueur de l'échangeur agit considérablement sur la résistance thermique. L'intensité
de la convection naturelle augmente en fonction du nombre de Rayleigh. Généralement,
on remarque que pour un Ra donné, plus la longueur du stratifié est importante, plus
ïmax est faible. Ce phénomène, similaire à celui observé par [16], est particulièrement
prononcé pour de faibles épaisseurs (L G [0.1 — 0.4]).

Après un certain point, l'augmentation de la longueur n'influence plus la résistance

thermique. En effet, la température maximale demeure fixe, malgré qu'il y ait davantage
de surface disponible pour l'échange thermofluide. En d'autres termes, le stratifié atteint
une longueur limite au-delà de laquelle toute matière additionnelle ne contribue plus
au transfert de chaleur. Dans ces conditions, une couche limite thermique se forme sur
la plaque chaude. Les longueurs efficaces observées correspondent donc à l'épaisseur de
la couche limite thermique à la sortie de l'échangeur, tel que :

Leff = Sr(y=l) (5.1)

où Leff désigne la longueur efficace. Numériquement cette limite est évaluée selon :

( Tmax [LeîI + A L j - T m a x ( L e / / j j / T m a x [I/j

AL -< o {.)

°ù Tmax (L) est la température maximale pour un L donné. Ainsi, lorsque l'amélioration
relative à l'augmentation de la longueur devient inférieure à 1%, on considère que la
longueur efficace est atteinte. Graphiquement, à la figure 5.1, cette limite, identifiée
par un carré noir, marque le début du plateau formé par la courbe Tmax(L). Pour des
14 18
valeurs de Ra variant de 10 à 10 , les longueurs efficaces, ainsi que les températures
maximales associées, sont présentées au tableau 5.2.

Les valeurs de Leff présentés au tableau précédent décroît avec l'augmentation du

 69

10

o
CD 10:

.iH L = 5.75
max
03
e oo e ao-0 »o£) e oo
0
3.5
IQDODD
-M 10
tri - Tni|1I'.||-|UJii^i u i i

-CD
-*- Ra = 1014
1.0 Ra = 1015
s - e
eu Ra = 1016
H 0.6 —*- Ra = 1017
Ra = 1018
10"
2 3 4 5 _ 6
L o n g u e u r du stratifié, L

FiG. 5.1 - Température maximale dans le stratifié en fonction du nombre de Rayleigh

et de la longueur,TV = 10, em = ec = 0, D = 0.001

Rayleigh. Ce phénomène indique que la couche limite thermique devient plus étroite
avec un Ra croissant. Cette observation est tout à fait cohérente avec les équations
d'approximation en régime couche limite (éq. 2.33-2.42). Selon l'approche analytique,
présentée à la section 2.3.2, l'épaisseur de la couche limite thermique à la sortie est :

ST = 1.65 Ra-Wk-WV'3

Soulignons que l'expression analytique précédente s'applique à une seule couche

uniforme opérant en régime couche limite thermique.

Lors d'une optimisation sans constrainte, l'AG converge vers une structure faite en­
tièrement de cuivre, à l'image du stratifié présenté au tableau 5.1. La performance ther­
mique devient alors uniquement dépendante de l'architecture poreuse. À la section 5.1.3,
la distribution de porosité au sein des structures optimales opérant en régime couche
limite thermique est abordée spécifiquement. Mentionnons que les stratifiés employés
pour construire la figure 5.1 comportent tous une distribution de porosité similaire.
 70

Globalement, on retrouve un milieu dense près de la paroi chaude et, en s'éloignant de

la source de chaleur, la porosité augmente graduellement.

Après avoir déterminé la longueur efficace du système comprenant 10 couches,

l'étape suivante consiste à déterminer si un refroidisseur non contraint, fabriqué à
partir d'une juxtaposition de plusieurs couches,- procure de meilleures performances
thermiques qu'une structure faite d'une seule strate uniforme. Ceci fait l'objet de la
prochaine section, où l'on s'intéresse à l'effet du nombre de couches du stratifié sur la
température maximale.
 71

TAB. 5.1 - Caractéristique et structure du design optimal pour : Ra = 1016,L = 0.5,N =

10, em = ec = 0, D = 0.001

Propriétés du stratifié
f M C
^ max
1.810E-4 1.282E+3 2.423E+4

Structure du stratifié
couches 1 2 3 4 5 6 7 8 9 10
matériaux eu eu eu eu eu eu eu eu eu eu
porosités 0.51 0.51 0.55 0.58 0.61 0.65 0.71 0.78 0.84 0.94

TAB. 5.2 - Longeur efficace en fonction du nombre de Rayleigh, N = 10, em = ec = 0,

D = 0.001
Ra 1014 1015 1016 1017 1018
L
eff 5.75 3.5 2.0 1.0 0.6
f1 7.56£ - 4 3.48£ - 4 1ME-4 7A7E - 5 3A6E - 5
- - max,num
 72

5.1.2 Effet du nombre de couches

Une caractéristique importante de l'architecture du stratifié est le nombre de couches

que compte ce dernier. Cette section porte donc sur l'impact du nombre de strates sur
la performance thermique d'un stratifié non contraint. Pour ce faire, plusieurs optimi­
sations sont réalisées en conservant les mêmes paramètres à l'exception du nombre de
sous divisions du système. En augmentant la complexité de la structure, on observe
une légère diminution de la température critique. Ce phénomène est perceptible pour
les différents Ra étudiés. L'influence de ce paramètre de design sur la distribition de
chaleur est nettement plus prononcée pour une variation entre N = 1 et N = 2. Géné­
ralement, on constate que l'influence du nombre de strates sur la température critique
s'amenuise rapidement avec l'augmentation de N, il s'agit donc d'une relation à ren­
dement décroissant. D'ailleurs, ces résultats concordent avec les travaux de [17], dans
lesquels on s'intéresse à l'effet du nombre de couches sur la résistance thermique d'un
empilement de micros conduits dans un contexte de refroidissement de composantes
électroniques. Les auteurs démontrent qu'il est profitable thermiquement d'accroître la
complexité de la structure, mais qu'après un certain point cela n'a plus d'impact sur le
transfert de chaleur.

À première vue les avantages que procure la stratification de la structure poreuse

semblent faibles. Cependant, dans certains contextes, ces gains peuvent tout de même
être recherchés. De plus, il faut préciser que ces résultats s'appliquent à une optimisation
sans contraintes. Le potentiel d'un échangeur fait de plusieurs strates se révèle lors de
la recherche d'un design optimal soumis à différentes contraintes. En effet, la flexibilité
de la constitution de l'échangeur, accordée par le nombre de couches, permet de mieux
s'adapter aux exigences de conception.

5.1.3 Distribution de porosité

Abordons maintenant la distribution de porosité au sein du stratifié à la suite d'une

optimisation sans contraintes. Globalement, on observe que le nombre de Rayleigh
affecte la répartition de matière au sein de l'échangeur. Pour de faibles nombres de
Rayleigh, l'algorithme génétique converge vers des architectures plus denses que pour
de grands Ra. Ce phénomène est illustré à la figure 5.2, où les résultats d'optimations
 73

réalisés pour L = 0.25 sont présentés. On souligne que d'après les longueurs effectives
présentées au tableau 5.2 un stratifié dont L = 0.25 n'opère pas en régime couche
limite, et ce, pour la gamme de Ra étudiée. Ainsi, lorsque le paramètre Ra est grand,
le fluide se dilate fortement et la vitesse dans les pores devient plus importante, ce qui
favorise le transfert thermique. En dimininuant le Rayleigh, l'échange de chaleur par
convection s'amenuise, alors le refroidisseur requiert une plus grande surface d'échange.
C'est ce qui explique, à la figure 5.2, la grande porosité du stratifié pour Ra = 1014
comparativement à la distribution de matière associée à Ra = 10 18 .

T r

0.9
"a
o
■©-
0 8
£ -
a
'â 0.7
o
*<D
4->
• rH
M
g 0.6
o
OH
0.5

0.4 J L

0 0.05 0.1 0.15 0.2 0.25

Position dans le stratifié, x

FlG. 5.2 - Distribution de porosité dans le stratifié en fonction de x, L = 0.25, N =

10,e m = ec = 0, D = 0.001

De plus, on remarque que la porosité augmente graduellement en s'éloignant de

la source de chaleur. Cette distribution de matière solide est analogue à celle décrite
dans les travaux de [49] qui portent sur le design d'une ailette en convection naturelle.
L'auteur démontre que pour des matériaux conducteurs, la forme idéale est de type
sharp-pointed. Cette géométrie comprend une concentration de matière solide à la base,
près de la source de chaleur. En s'éloignant de la base, ce design optimal offre une plus
grande portion d'espace destinée au fluide refroidissant. Similairement, l'architecture
poreuse du stratifié, obtenue à la suite de l'optimisation par algorithme génétique,
comprend des couches poreuses qui se densifient à l'approche de la paroi chaude. Cette
 74

configuration répartit adéquatement la chaleur au sien de l'échangeur, ce qui permet

d'exploiter le plein potentiel de l'espace d'échange disponible.

On a fait une autre série de simulations avec de longs échangeurs afin de vérifier
la distribution de porosité quand le système opère en régime couche limite thermique.
Pour ce type d'écoulement, on observe que la distribution de porosité évolue de manière
similaire peu importe le nombre de Rayleigh considéré. La figure 5.3 illustre la porosité
des couches en fonction de la position relative dans la couche limite. Les valeurs de
longueurs effectives, obtenues à la section 5.1.1 sont employées en guise de valeurs de
références. L'abscisse de cette figure se définie comme : x/LefftRa. Par exemple, pour
lire une valeur correspondante à x — ÔT, il faut se référé à l'abscisse 1.0. Grâce à
cette représentation, il est possible de comparer sur une base commune l'architecture
poreuse des stratifiés soumis à différents Ra. D'après ce graphique, on remarque une
certaine ressemblance entre les optimisations. En effet, malgré que l'épaisseur ST varie
en fonction de l'intensité du Ra, les distributions de porosités au sein de la couche
limite sont semblables. Pour les couches adjacentes à la plaque chaude, l'optimisation
prescrit une porosité de 0.414. Ce résultat concorde avec les valeurs escomptées d'après
l'approche analytique présentée à la section 2.3.3. Dans la deuxième moitié du stratifié,
la densité de matière solide diminue progressivement pour atteinte de grandes porosités
à x = LefftRa.
 0.9
-*- Ra = 1014 *
1
0.85 - G- Ra = 1015 O
D Ra = 1016 //
0.8 —*— Ra = 1017 l'f -
0.75 Ra = 1018 i 7/
0
o

O
o!'
CS
e
• r - l
4->
0.7
0.65
0.6 / */ ; ri
/
/Pi à"
-0) / / 0 .' / f
-M
0.55 - ' /
• r-H
^ / ■ D

/ /

c .#
O 0.5 -* /
PLH -0 ' ,. J
s- *'
0.45 #— ■ —Mf- rO^
^ --•
— .^
0.4 1 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Position d a n s le stratifié, x

FlG. 5.3 - Distribution de porosité dans le stratifié en fonction de x/Leff^a, N =

10,e m = ec = 0, D = 0.001

5.1.4 Masse et coût

L'optimisation sans contraintes, on le rapelle, converge invénitablement vers un

échangeur constitué entièrement de cuivre, car, bien que ce matériau soit le plus cher et
le plus massif parmi les constituants possibles (tab. 2.1), il détient la meilleure conduc-
tivité thermique. Dans ce cas, la masse et le coût des designs optimum dépendent
entièrement de la distribution de porosité telle que présentée à la section précédente
(sec. 5.1.3). De plus, étant donné que l'ensemble des designs sont faits du même ma­
tériau, la variation du coût d'un échangeur en fonction de la longueur correspond à
la variation de la masse (figure 5.4) mulpliée par le prix adimensionnel du cuivre. Sur
cette figure, on remarque que la masse du stratifié croît rapidement avec l'augmentation
de longueur pour des valeurs de L comprises entre 0 et L e / / / 2 . C'est d'ailleurs dans
cette plage de longueurs que la température critique subit les plus fortes diminutions.
Pour des dimensions supérieures à la moitié de la longueur efficace, l'accroissement de
la masse est moins prononcé, et la baisse de T max est aussi moins importante.

Il est possible de comprendre la variation de la masse globale du stratifié grâce aux

 70

10 6

■ ; ;

( >
10 4 - * " *
o - -O" nu • ■ ■ ■ 1.1
» -Cr „.^ .n- _n, ■

10 3

o
CD
M
<s i o 2 -#-■ Ra=1014
- 9 - Ra = 1015
D Ra=1016
—*— Ra = 1017
-- - Ra = 1018
10 1
0 0.2 0.4 0.6 0.8
L o n g u e u r du stratifié, L / Leff, Ra

FlG. 5.4 - Masse du stratifié en fonction de la longueur et du nombre de Rayleigh,

AT = 10, em = ec = 0, D = 0.001

distributions de porosité associées aux différentes longueurs d'empilement. Par exemple,

la figure 5.5 présente l'arrangement poreux des stratifiés optimums obtenus avec un
nombre de Rayleigh de 1016 pour différentes longueurs de L. D'après ce graphique, on
constate que l'échangeur se densifie progressivement en augmentant sa longueur. Avec
l'accroissement de la taille du refroidisseur, les porosités des couches adjacentes à la
source paroi chaude tendent vers la porosité optimale </>opi = 0.414 afin de favoriser la
répartition de chaleur au sein de la structure. Ceci a pour effet d'accroître le poids de
l'empilement.
 77

1 1 i i i i i i

f ■-•*-• L = 0.25

* / p - e - L = 0.50
0.9
' / n L = 0.75
i é / '
□ / 1 —*—L=1.00
tf 0-8
_ é 1 t -
F 1 f é / ' 1
eg L=1.50
f !
a * 0 D / ' ' ■-*-■ L = 2.00
U 0.7 ; / 7
'

*->
o 7 / / /
/ / 7
'§ 0.6 * ©-0 P / /'
O
fin i ' y
/ ' /
0.5 i D
0-0 .□ j /
é a ET - ^* *^
0- ■□' ET
■ — * ■ "

0.4 1 1 1 1 1 1

0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.

Position dans le stratifié, x

FlG. 5.5 - Porosité en fonction de la longueur, Ra — 10 ,e m = ec — 0, D = 0.001

5.2 Optimisation sous une contrainte de niasse

Après avoir étudié les résultats de l'optimisation sans contraintes du stratifié poreux,
l'étape suivante consiste à imposer une contrainte au système. Dans un premier temps,
on observe l'effet d'une masse maximale permise Mco sur la recherche d'un design
optimal. La valeur de Mco est déterminée à partir des résultats de l'optimisation non
contrainte obtenue à la section précédente, car pour que cette limite soit effective, il
faut que la masse maximale imposée soit inférieure à celle obtenue dans le cas non
contraint. Ainsi, les limites de masse correspondent à une certaine fraction des valeurs
présentées à la figure 5.4. Lors des simulations, le même nombre de variables de design
est octroyé au système. L'échangeur comprend 10 couches qui peuvent être fabriquées
à partir des matériaux énumérés au tableau 2.1. La précision des variables reliées à
la porosité demeure inchangée (voir section 4.2.2). On rappelle que l'objectif est de
minimiser la température maximale dans le stratifié, T m a x , tout en tenant compte d'une
masse limite. Dans ce cas, la fonction objectif correspond à l'équation (4.2) qui après
 78

simplification s'écrit :

M
/ (x) = fmax \l + eT max -1,0
Mm

L'optimisation est réalisée pour un nombre de Rayleigh de 10 15 . La longueur du

stratifié étudié est L = 1.0, ce qui correspond à la moitié de la longueur effective,
comme l'indique le tableau 5.2. La figure 5.6 illustre l'architecture poreuse ainsi que les
principales caractéristiques du stratifié optimum obtenu pour le cas non contraint. Ces
résultats serviront de référence afin de mesurer l'impact de la contrainte de masse sur la
structure optimale. Les résultats de simulations présentés à la figure 5.7 sont réalisées
avec des masses critiques Mco égale à 75%, 50%, 25% et 10% de la valeur de la masse
calculée lors l'optimisation sans contraintes (fig. 5.6). Les résultats sont présentés à la
figure 5.7.

M c o /M s c o =0,C c o /C 8 C O =0
i
o/)
0.8
o./
0.6
o.',
0.4
0.3 T = 3.960E -4
I 0.2
0.1
max
M =2.488E+3
C = 4.703E+4
1 2 3 4 5 6 7 8 9 10
Couches du stratifié

F I G . 5.6 - Structure poreuse optimale sans contrainte, Ra — 10 15 ,L 1.0,N = 10,

M/Msœ = 0, C/Csco = 0, D = 0.001

Dans chacun des cas étudiés, la contrainte de masse est respectée. En fait, on
constate que la masse des designs présentés à la figure 5.7 correspond à la contrainte
imposée. Dans l'espace de design, les structures optimales se situent donc à la fron­
tière de la contrainte. Afin de réduire le poids du stratifié, tout en mininimisant la
température maximale, l'AG emploie trois mécanismes : la substitution de matériaux,
l'augmentation de porosité et l'élimination de couches.

Pour une contrainte de masse de 75% (optimisation 1, fig 5.7), on observe que l'AG
remplace certaines couches de cuivre par de l'aluminium, un matériau qui offre une
bonne conductivité thermique tout en étant relativement léger. Les couches modifiées
 79

Optimisation 1 Optimisation 2
M c o / M s c o = 0.75, e c =0 M
co /M
sco=0-50>Ec:
i 1
0.9 h 0.9
0.8 0.8
-eu 0.7 0.7
+J
'S 0.6 0.6
o
S 0.5 0.5

0.4 0.4

0.3 T = 4.008E -4 0.3 T =4.162E-4

0.2 max
0.2 M = 1.870E+3
0.1 M = 1.252E+3
0.1 C = 3.469E+4
C = 2.162E+4
3 4 5 6 7 8 10 3 4 5 6 7 8 10
Couches du stratifié Couches du stratifié

Optimisation 3 Optimisation 4
M C O /M S C O =0.25, e c =0 Mco/MSco=0'10.Eo=0
i 1
0.9 0.9
0.8 0.8
0.7 -fl> 0.7
+J

0.6 0.6
o
0.5 '-A
0.5
r,
0H
0.4 0.4
0.3 T = 4.505E-4 0.3 T =5.185E-4
max max
0.2 0.2
M =6.219E+2 M =2.537E+2
0.1 0.1
C = 9.286E+3 C = 3.910E+3
1 2 3 4 5 6 7 8 9 10 2 3 4 5 6 7 8 9 10
Couches du stratifié Couches du stratifié

FiG. 5.7 - Structures poreuses optimales sous une contrainte de masse, Ra = 1015,
L - 1.0, TV - 10, Mœ/Msœ = [0.75, 0.5, 0.25, 0.1], Cco/Csco = 0, D = 0.001
 80

se situent à l'opposé de la plaque chaude. Comme il a été démontré à la section 5.1, les
échanges thermiques sont de forte intensité près de la plaque chaude et s'amenuisent
en s'éloignant de la source de chaleur. Par conséquent, les premières strates sont de
grande importance au plan du transfert thermique et doivent être conçues de manière à
distribuer efficacement la chaleur au sein de l'échangeur. Pour ce faire, les couches adja­
centes à la source de chaleur requièrent un matériau possédant une grande conductivité
thermique. Ainsi, le fait de cibler les dernières couches pour réduire la masse globale
de l'échangeur s'avère être un choix judicieux. On observe à la figure 5.6 que les résul­
tats issus de l'AG vont en ce sens. La première portion du stratifié, celle qui participe
le plus activement au transfert de chaleur, est faite de cuivre, ce matériau présentant
la meilleure conductivité thermique parmi les choix disponibles. La combinaison al-cu
offre un bon compromis entre la masse de la structure et les propriétés thermiques du
refroidisseur. D'ailleurs, on retrouve ce type de combinaison dans plusieurs modèles
d'échangeur de chaleur dédiés au refroidissement de composants électroniques [50].

En accentuant la contrainte de masse à 50%, l'ajout d'aluminium dans la structure

devient plus important. Pour l'optimisation 2, l'échangeur contient autant de couches
de cuivre que de couches d'aluminium. On note que les couches adjacentes à la source
de chaleur sont toujours faites de cuivre afin de favoriser le transfert de chaleur. De
plus, on remarque que les porosités sont globalement plus élevées que dans le cas sans
contraintes présentées à la figure 5.6. Cette modification réduit la quantité de matière
solide, ce qui abaisse la masse du stratifié. De cette façon, une portion de l'échangeur
peut être faite de cuivre tout en respectant la masse maximale.

Lorsque la limite massique devient encore plus importante, l'AG accentue davantage
la porosité globale du stratifié puis élimine des couches. Pour une contrainte de masse
de Mco/Msco =25%, la structure optimale est entièrement constituée d'aluminium, à
l'exception de la première couche qui demeure en cuivre. Afin de diminuer la masse de
l'échangeur, la porosité des couches peut être augmentée. Cependant, l'augmentation
de la proportion de fluide dans les strates détériore la conductivité équivalente keq,
telle que décrite à l'équation 2.14, et cela a pour effet de limiter le transfert de chaleur
par conduction au profit de la convection dans les premières couches. De manière à
conserver un juste milieu entre ces deux modes d'échange thermique, l'AG procède à
l'élimination des dernières couches. À la figure 5.7, les deux configurations qui résultent
de l'optimisation sous des contraintes de 25% et 10% de masse comprennent des couches
vides. Grâce à ce mécanisme, la longueur du stratifié est indirectement optimisée.
 81

L'ensemble des modifications réalisées pour respecter la contrainte imposée altère

les performances de refroidissement. En effet, l'augmentation de la valeur de T m a x pour
les contraintes de 75%, 50%, 25% et 10% est respectivement de 1.2%, 5.1%, 13.8% et
28.1% par rapport à l'optisation sans contrainte (fig. 5.6).

La répétitivité des résultats présentés à la figure 5.6 est particulièrement bonne.

Lors de cette étude, pour chacune des optimisations effectuées, cinq simulations avec des
paramètres identiques sont lancées, et les variations sur la valeur de Tmax sont inférieures
à 0.6%. Les différentes configurations produites représentent alors des systèmes quasi
optimaux. Pour chacun des cas, le meilleur des résultats est présenté à la figure 5.7.

5.3 Optimisation sous une contrainte de coût

À cette étape, aucune restriction n'est imposée à la masse du stratifié. Une contrainte
de coût Cco est plutôt appliquée au système. La procédure utilisée est similaire à celle
présentée à la section 5.2. Le coût critique correspond à une fraction de la valeur calculée
lors de l'optimisation sans contraintes (voir fig. 5.6). Comme ce fut le cas dans la section
précédente, les valeurs de 75%, 50%, 25% et 10% du coût d'une structure optimale non
contrainte sont sélectionnés. On emploie un stratifié identique à celui de l'optimisation
sous une contrainte de masse soit : Ra = 10 15 , L = 1.0 et N = 10. Les nouveaux designs
obtenus d'après la fonction objectif décrite à l'équation (4.2) et réécrite sous la forme
simplifiée suivante sont présentés à la figure 5.8.

/ (x) = r m a x | f ec max ( -~ 1,0

Dans le but de réduire le prix de la structure, l'AG emploie les mêmes stratégies
que lors de l'optimisation de la masse. Par exemple, dans le cas de la contrainte de
Cœ/Csco =75%, l'algorithme substitue deux des dernières couches de cuivre par de
l'aluminium, qui représente un matériau moins coûteux et qui possède de bonnes pro­
priétés thermiques. Le stratifié conserve donc une large portion de cuivre près de la
paroi chaude ce qui favorise la distribution de chaleur. La configuration de la structure
poreuse ainsi que les performances thermiques du stratifié issu d'une optimisation sous
 82

Optimisation 1 Optimisation 2
em=0,Cco/Csco=0.75 <V=0,C C O /C S C O =0.50
i
0.9
0.8
St
0.7
^ 0.6
g 0.5 3
o
I 0.4 h o 0.4

T = 4.005E -4 0.3 T =4.127E-4

max max
0.2
M = 1.898E+3 0.1 M = 1.398E+2
C = 3.521E+4 C =2.352E+3
3 4 5 6 7 8 4 5 6 7 8 9 10
Couches du stratifié Couches du stratifié

Op timi satil m 3 Optimisation 4

em=0,Cco/Csco=0.25 E m =0, C c o / C s c o = 0 . 1 0
1
0.9 St
St 0.8 Al St
Al Al St 0.97
- Al St 0.90 0.7 Al 0.90 0.90
O.IM 0.84 0.84
0.77 0.77 Al 0.77
Cu Al Al .•S
m
0.6 Al 0.71
Al Al Al

l
0.67 0.65 0.68 0.65
Al 0.58
0.61 S 0.5 0.61
0.55
° 0.4
ÇL, „.
f =4.374E-4 0.3 f =4.791E-4
max max
0.2
M = 1.005E+3 0.1 M = 5.528E+2
C = 1.177E+4 C = 4.715E+3
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10
Cou ches du strat ifié Couche s dvi stra tifié

FiG. 5.8 - Structures poreuse optimales sous une contrainte de masse, Ra L015
L = 1.0, A^ = 10, M/Msco = [0.75,0.5,0.25,0.1], C/Csco = 0, D = 0.001
 <s:i

une contrainte de coût de 75% sont similaires aux résultats obtenus à la suite d'une
optimisation sur la masse de Mco/Msco =75%. Ces deux approches conduisent donc à
la même région de l'espace de design.

Lorsque la contrainte de coût s'accentue, l'acier devient un matériau envisageable

en raison de son très faible coût. Ainsi, l'échangeur optimal doté de la moitié du coût du
modèle non contraint possède une couche d'acier positionnée à l'opposée de la plaque
chaude. Etant donné qu'il n'y a pas de contrainte de masse, l'AG préfère ajouter une
couche d'un matériau peu cher et peu performant thermiquement plutôt que d'éliminer
une partie du stratifié. De plus, le fait de choisir l'acier comme substituant permet
de conserver des couches de cuivre relativement denses près de la source de chaleur,
une région jugée de grande importance sur le plan du transfert thermique. Pour des
restrictions importantes sur le prix, comme dans le cas de la simulation réalisée avec un
coût critique de Cco/Csco = 10% , ce compromis ne fonctionne plus. L'algorithme doit
éliminer les dernières strates afin d'attribuer des matériaux et des porosités adéquats
aux premières couches. Aussi, on note que la porosité globale des stratifiés optimaux
issus d'une optimisation sous une contrainte de coût est légèrement plus faible que
lors de l'optimisation avec des contraintes sur la masse. Sur le plan de la température
maximale, les augmentations engendrées lors des d'optimisation sont inférieures à celles
produites pour des contraintes de masse. Les augmentations de T m a x par rapport au
cas non contraint sont de 1.1%, 4.1%, 10.3% et 21% pour les quatre cas présentés à la
figure 5.8.

5.4 Optimisation sous une combinaison de

contraintes de masse et de coût

Cette section se consacre à l'étude de l'influence d'une combinaison de contraintes

de masse et de coût sur l'architecture optimale du stratifié. De manière à bien mesurer
l'effet de cette double contrainte, on traite le même stratifié que celui utilisé aux sections
précédentes (Ra = 10 15 , L = 1.0 et N = 10).

Le prix du refroidisseur poreux est intimement lié à la masse, comme l'indique

l'équation (2.53). Dans le but de mieux comprendre l'interaction entre la contrainte de
 84

masse et celle de coût, plusieurs optimisations sont réalisées pour différents rapports
de Mco/Msco et Cco/Csco. Les résultats sont présentés à la figure 5.9. Également, des
exemples détaillés d'optimisations de combinaisons de contraintes (marquées par des
points jaunes) sont illustrés à la figure 5.10. Ces cas montrent l'allure générale des
stratifiés issus des principales zone d'optimisation que l'on retrouve à la figure 5.9.

o
o
IO Contrainte de
c
o
masse dominante
t;
O 0.6 — a — contrainte de masse
iO contrainte de coût
□ cas typiques
<3 O double contraintes
o
o
CD
-S
CD Contrainte de
coût dominante
-t-J
PI
o
o 5.530E A
J I
0.1 0.2 0.3 0.4 0.5 0.6

Contrainte de masse, M c o / M g c o

FlG. 5.9 - Optimisation de la structure poreuse sous une contrainte de masse et de

coût, Ra = IO 15 , L = 1.0, N = 10, D = 0.001

La première courbe, identifiée par des carrés bleus, correspond aux optimisations
exécutées seulement sous une contrainte de masse. On remarque que selon ces condi­
tions, la relation entre le poids et le prix de l'échangeur est pratiquement linéaire. La
deuxième courbe, tracée en rouge, concerne les résultats issus d'optimisations impli­
quant uniquement une contrainte de coût. Comme pour la première approche, la masse
du stratifié varie proportionnellement au prix imposé lorsque les contraintes de coût
Cco/Csœ sont supérieures à 50%. Dans cette région, les deux courbes sont d'ailleurs très
près l'une de l'autre, ce qui indique qu'il est possible d'obtenir un même design par l'en-
 85

tremise des deux types de contraintes. En d'autres termes, si on réalise une optimisation
(dans la région où les deux contraintes se superposent) sous une contrainte de masse, la
structure optimale aura un coût associé. Si on utilise ce coût comme contrainte unique,
on obtiendra de nouveau le même design. En diminuant le coût maximal admissible, la
différence entre les configurations optimales issues des contraintes de masse et de coût
s'accentue. Les optimisations 3 et 4 de la figure 5.10 sont un bel exemple de cette diver­
gence. Ce phénomène s'explique par l'apparition de couches d'acier dans la structure.
Pour un prix donné (Cco/Csco), la sélection de l'acier dans la fabrication produit des
stratifiés beaucoup plus massiques comme le témoigne l'optimisation 4 de la figure 5.10.
Effectivement, dans le cas d'une contrainte de masse, les designs optimaux sont tous
constitués de cuivre et/ou d'aluminium (voir optimisation 1, 2 et 3 à la figure 5.10), car
les strates d'acier alourdissent grandement la structure compte tenu du gain thermique
qu'elles procurent. Lorsque l'accent est plus porté sur le faible coût de l'échangeur, les
stratifiés optimaux sont généralement plus longs et plus denses. En somme, la figure 5.9
met en relief trois régions distinctives. Dans le cas où l'on réaliserait une optimisation
avec simultanément une contrainte de masse et une de coût, si on se situe au dessus de
la ligne bleu, seule la contrainte de masse sera effective. De même, sous la ligne rouge,
seule la contrainte de coût est active. Finalement dans la zone entre les deux lignes, on
retrouve les combinaisons de contraintes menant à deux contraintes actives.

Afin de mettre en évidence des cas où les deux contraintes seraient actives simul­
tanément, les limites de masse et de coût sont sélectionnées judicieusement. Quelques
optimisations sont réalisées à l'intérieur de la région formée par les courbes de contrainte
de masse et de coût. Les combinaisons simultanées de contraintes identifiées par des
cercles verts sur la figure 5.9 et détaillés à la figure 5.11 indiquent qu'il est possible
d'imposer deux types de restrictions effectives. En effet, les valeurs prescrites de masse
et du coût sont plafonnées simultanément. Afin de réduire le poids global, les échan-
geurs de la figure 5.11 sont majoritairement fabriqués d'aluminium. On remarque que
pour répondre aux exigences demandées l'AG emploie deux types de combinaisons de
matériaux. L'association aluminium-acier réduit le coût tandis que la combinaison al-
cu minimise plutôt la masse. Fait intéressant à noter, le stratifié issu des contraintes
Mco/Msco = 30% et Cco/Csco = 20% est constitué des trois matériaux disponibles. De
plus, la position des matériaux dans ce refroidisseur concorde avec les observations réa­
lisées aux sections précédentes (section 5.2-5.3). Le cuivre se situe près de la source de
chaleur alors que l'acier est relégué à la dernière portion de l'échangeur.
 86

Optimisation 1 Optimisation 2
M c o / M s c o = 0.95, C œ / C s c o = 0.95 M n n / M c „ =0.50, C _ / C =0.50
i 1
0.9 0.9 ■
0.8 0.8
0.7 0.7
-V 0.6 - £ 0.6

o 0 5 '% 0.5 h
s*
o 0.4 O 0.4
0.3 T = 3.961E -4 0.3 T = 4.162E -4
0.2 max 0.2 max

0.1
M =2.362E+3 0.1 M = 1.246E+3
C = 4.464E+4
C = 2.091E+4
2 3 4 5 6 7 10 2 3 4 5 6 7 10
Couches du stratifié Couches du stratifié

Optimisation 3 Optimisation 4
M c o / M s c o = 0.20, C c o / C g c o = 0.90 M c o / M s c o = 0.90, C c o / C s œ = 0.20
i i
0.9 0.9
0.8 Al 0.8
0.7 Al OIM 0.7
^ 0.6 Al (17-1 •g 0.6
| 0.5 0.65 '§ 0.5
SH
O 0.4 O 0.4
CM P-l
0.3
T = 4.629E -4 0.3
T = 4.474E -4
max max
0.2 0.2
M = 4.967E+3 0.1
M =7.255E+3
0.1
C = 6.358E+4 C = 9.451E+4
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10
Couches du stratifié Couches du stratifié

FiG. 5.10 - Optimisation de la structure poreuse sous une combinaison de contraintes

de masse et de coût, Ra = 1015, L = 1.0, iV = 10, M/Msco = [0.95,0.5,0.2,0.9],
C/Csco = [0.95, 0.5,0.9, 0.2], D = 0.001
 87

Optimisation 1 Optimisation 2
M
co /M
s c o = 0.20, C c o / C s c o = 0.1 1
M™, / MQ„ = 0.30, C„„ / C c = 0.2
i
0.9 ■ St 0.9
0.8 0.8

-V 0.7 MB 0.7
-M
' S 0.6 -g 0-6
g | 0.5
g 0.5 h
PH 0.4

T = 4.797E -4 0.3 T = 4.486E -4

max
max 02
M =4.715E+2 0.1
M = 7.489E+2
C = 4.724E+3 C = 9.356E+3
3 4 5 6 7 10
3 4 5 6 7 10

Couches du stratifié Couches du stratifié

Optimisation 3 Optimisation 4
M c o / M s c o = 0.35, C c o / C s c o = 0.15 , M c o / M s c o = 0.35, C c o / C s c o = 0.3
0.9 Al
St Al
0.8 0.97
Al St 0.94
Al 0.90
Al .<U 0.7 ÛJ Al Al 0.84
Al 0.84 St 0.84 -M Ol Al
• 0.77 3 0.6 0.74 U./4
Al AI 0.74 0.74 0.70
Al Al 0.68
Al 0.65 og 0.5 0.61
Al 0.58 0.61 0.58
0.55 ^ 0.4
f =4.568E-4 0.3 f =4.321E-4
max
0.2 M = 8.690E+2
M =8.559E+2
0.1
C = 7.065E+3 C = 1.368E+3
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10

Coucr tes d U S t ]•atifi é Ce>uch 3S d li str atifié

FiG. 5.11 - Optimisation de la structure poreuse sous une combinaison de contraintes

de masse et de coût, Ra = 10 15 , L = 1.0, N = 10, M/Msco = [0.2,0.3,0.35,0.35],
C/Csco = [0.1,0.2,0.15,0.3], D = 0.001
 À la lumière de ces résultats, on constate que la constitution d'un stratifié opti­
mal est fortement dictée par la paire de contraintes de masse et de coût. On remarque
d'ailleurs certaines tendances quant à la composition des échangeurs selon les régions
du domaine formées par les contraintes Mco/Msco et Cco/Csco. Afin de mieux perce­
voir l'arrangement optimal en fonction des limitations imposées, des optimisations sont
réalisés sur l'ensemble du domaine. Les figures 5.12, 5.13 et 5.14 présentent respective­
ment le nombre de couches de cuivre, d'aluminium et d'acier dans la structure poreuse
optimale selon les combinaisons de contraintes.

FlG. 5.12 - Nombre de couches de cuivre dans la structure poreuse optimale selon la
combinaison de contraintes de masse et de côut, Ra = 1015, L = 1.0, N = 10, D — 0.001

Aux figures 5.12 à 5.14, le gradient de couleur indique la teneur du matériau sélec­
tionné au sein de l'échangeur. Par exemple, les designs constitués majoritairement de
cuivre se concentrent dans la région où les contraintes sont moins importantes (coin
supérieur droit). En accentuant les contraintes, le nombre de couches de cuivre di­
minue progressivement. Près de l'origine, endroit du graphique où la combinaison de
contraintes est très restrictive en terme de masse et de coût, le cuivre est complètement
exclu de la composition des stratifiés optimaux. En ce qui concerne l'acier (figure 5.14),
 89

!U

lu

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

M C o/M s c o

FiG. 5.13 - Nombre de couches d'aluminium dans la structure poreuse optimale selon
la combinaison de contraintes de masse et de côut, Ra = 1015, L = 1.0, N = 10,
D = 0.001
 !)()

!U

lu

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

M c o /M s c o

FiG. 5.14 - Nombre de couches d'acier dans la structure poreuse selon la combinaison
de contraintes de masse et de côut, Ra = 1015, L = 1.0, N = 10, D = 0.001
 !)l

i r~

!U

!U

_i i ' i_

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

M c o /M s c o

FiG. 5.15 - Nombre de couches vide dans la structure poreuse selon la combinaison de
contraintes de masse et de côut, Ra = 1015, L — 1.0, N = 10, D = 0.001
 92

son application se limite à la portion du domaine située sous la courbe des contraintes de
masse présentée à la figure 5.9. La plus forte utilisation de l'acier dans la structure opti­
male se produit sous d'importantes contraintes de coût jumelées à de faibles limitations
sur la masse du système. L'aluminium pour sa part est fréquemment employé comme
le témoigne la figure 5.13. Ce matériau est d'ailleurs employé à la fois pour diminuer
le coût et réduire la masse. On retrouve donc des couches d'aluminium dans plupart
des architectures optimales issues de problèmes contraints. L'aluminium est combiné
au cuivre pour abaisser la masse de la structure et est jumelé à l'acier pour atteindre de
faibles coûts. On remarque aussi que l'emploi de l'alumium est plus prononcé lorsque les
limitations imposées sont importantes. Ce phénomène est illustré à la figure 5.13 par les
régions foncées. Lorsque que les limitations deviennent trop importantes, l'AG n'a pas
d'autre choix que de tronquer le stratifié. Dans ce cas, aucun matériau n'est attribué
à une ou plusieurs couches de l'échangeur. Comme il a été démontré à la section 5.2,
les couches vides se situent à l'opposé de la source de chaleur. La figure 5.15 met en
relief les régions du domaine de contraintes où la longueur de l'échangeur est réduite.
Dans le cas étudié {Ra = 10 15 , L = 1.0), l'élimination de strates se produit pour des
contraintes de masse de l'ordre de Mco/Msco < 40%. En somme, les figures 5.13 à
5.15 permettent de visualiser rapidement la composition générale d'un stratifié selon
les contraintes sélectionnées. Ces graphiques, de pair avec les distributions optimales
de porosité présentée à la section 5.1.3, forment un portrait complet de l'architecture
optimale de l'échangeur poreux sous à une combinaison de contraintes de masse et de
coût.

L'addition d'une deuxième contrainte conduit à des solutions quasi identiques à

l'intérieur de la région délimitée par les courbes de contrainte de masse et de coût.
Globalement, au sein de cette zone, la température maximale augmente graduellement
en se déplaçant vers l'origine. Cette conclusion est cohérente avec le fait que le domaine
de design se réduit lorsque les contraintes s'accentuent. Au plan de la répétabilité des
résultats, la variation des valeurs de Tmax obtenues pour plusieurs simulations identiques
demeure toujours inférieure à 1% dans l'ensemble des essais réalisés ici.
 Conclusion

Dans cette étude, l'optimisation d'un stratifié poreux opérant en convection natu­
relle a été réalisée. Les résultats indiquent qu'il est possible d'améliorer les performances
thermiques du refroidisseur en sélectionnant une distribution de porosités adéquate ainsi
qu'un ordonnancement idéal des matériaux de la série de couches. Les travaux révèlent
également que l'on peut déterminer des structures présentant de bonnes propriétés ther­
miques tout en respectant des contraintes reliées à la masse et au coût. Un des aspects
novateurs de cette recherche est l'emploi d'un algorithme génétique dans le processus
d'optimisation. Cette technique permet de cerner rapidement un optimum d'un vaste
domaine en évaluant seulement une faible portion des solutions potentielles. Afin de
coupler le modèle numérique à l'AG, le problème fort complexe que représente l'optimi­
sation du système thermofluide a été ramené à une expression simple qui tient compte
de la température, de la masse et du coût. Cette approche permet d'établir une in­
teraction entre l'objectif de la recherche et l'architecture du refroidisseur. La méthode
employée dans ce document se distingue des autres travaux réalisés sur le sujet quant
à l'optimisation simultanée du type de matériaux et de la porosité.

Les résultats concluants obtenus à la suite de la présente étude ouvrent la porte

à des recherches plus approfondies sur le sujet. Il serait intéressant de considérer un
modèle mathématique plus complexe tenant compte de la conduction axiale et du rayon­
nement, par exemple grâce à l'approximation de Rosseland [51]. En effet, à de fortes
températures, le rayonnement thermique devient important et affecte significativement

93
 94

le transfert de chaleur et la distribution de température [52]. De plus, le cas d'un milieu

poreux sans équilibre thermodynamique local pourrait être traité [53, 54]. Cette modifi­
cation permettrait de considérer un modèle où, en un point du domaine, la température
du fluide et celle de la matrice solide ne sont pas équivalentes. Ceci impliquerait qu'il
faille résoudre l'équation d'énergie pour la partie solide et la partie fluide indépendam­
ment, tout en tenant compte de leur interaction réciproque. Aussi, dans ce mémoire,
le système étudié opère en régime permanent. Il serait pertinent d'envisager l'aspect
temporel dans la résolution thermofluide car, du point vue pratique, on retrouve plu­
sieurs applications industrielles impliquant une source de chaleur variable [55]. En ce qui
concerne l'écoulement de fluide, un modèle plus complexe qui considère les effets reliés
à la turbulence mériterait d'être envisagé. Afin de tenir compte des pertes de charge
associées à la turbulence du fluide dans les pores, le modèle de Forcheimer [34, 56]
pourrait être greffe à la loi de Darcy. Cette dernière modification, jumelée à l'addition
d'un modèle radiatif, permettrait de considérer des stratifiés poreux opérants à grand
Rayleigh.

Sur le plan de l'architecture interne du milieu, il serait intéressant d'examiner l'effet

de la géométrie des micro-conduits sur le transfert thermique. La forme des canaux pour­
rait même être optimisée. En ce sens, plusieurs structures poreuses présentant divers
types de perméabilités pourraient être ajoutées à la banque de matériaux disponibles
lors de l'optimisation. Alors, en plus de déterminer la porosité idéale, l'AG choisirait
le type de pores qui favorisent la dissipation de chaleur, tout en limitant la perte de
charge. Au cours du présent travail, le stratifié étudié possède une forme rectangulaire
pouvant varié de longueur grâce au rapport L. Dans l'optique d'optimiser la masse, le
coût et même l'espace, il serait à propos de considérer d'autres formes d'échangeurs.
Cet exercice pourrait d'ailleurs conduire à une optimisation conjointe de la géométrie
externe et interne du refroidisseur.

Relativement à l'aspect pratique de la fabrication d'un échangeur poreux, l'indice de

performance d'un système devrait tenir compte des contraintes reliées à la production
du design optimal. Par exemple, le coût d'un modèle de refroidisseur serait fonction à la
fois du prix des matériaux sélectionnés et des coûts associés à l'assemblage du système.
Dans certains appareils, on retrouve plusieurs échangeurs de chaleur, au lieu d'optimiser
chaque refroidisseur individuellement, il serait pertinent de considérer l'ensemble du
système. Dans ce contexte, l'optimisation du contrôle thermique aborderait globalement
des aspects comme la maintenance et les coûts de production.
 95

Toutes ces propositions complexifleraient le modèle et, du même coup, élargiront

l'espace de design. Par conséquent, les algorithmes génétiques employés au cours de ce
mémoire devront possiblement être améliorés afin d'identifier efficacement les solutions
optimales. Pour ce faire, un AG travaillant avec de multiples populations pourrait être
implémenté (rnulti-population genetic algoritms MPGA). Cette stratégie consiste à
fractionner les solutions potentielles en sous-groupe qui évoluent indépendants et dans
lesquels les opérateurs génétiques conventionnels sont appliqués. Ainsi, chaque popula­
tion est libre d'emprunter un parcours différent dans l'espace de design. À l'occasion,
des échanges génétiques peuvent être possibles entre les groupes grâce à la migration
de certains individus, permettant ainsi de mieux orienter la recherche d'un optimum.
Cette approche améliore l'efficacité de la recherche en maintenant une diversité globale
à travers les populations et prévient la convergence prématurée [57, 58].

En conclusion, la convection naturelle en milieu poreux représente un champ de

recherche prometteur en raison de ces multiples applications potentielles en ingénierie
telles que : le refroidissement de composante électronique, les collecteurs solaires, l'isola­
tion thermique, la climatisation, la ventilation, le refroidissement de réacteur nucléaire
et la géothermie [34].
 Bibliographie

[1] G. E. Moore, Cramming more components onto integrated circuits, Electronics 38

(8)(1965).
[2] B. Yu, M. Meyyappan, Nanotechnology : Rôle in emerging nanoelectronics, Solid-
State Electronics 50 (4) (2006), 536-544.
[3] D. S. Steinberg, Cooling Techniques for Electronics Equipment, 2nd Ed, John
Wiley and Son inc, 1991.
[4] A. McWilliams, The market for thermal management thechnologies, Business Com­
munications Company, GB-SMC024D (2006), 281.
[5] F. Bobaru and S. Rachakonda, Optimal shape profiles for cooling fins of high and
low conductivity, International Journal of Heat and Mass Transfer 47 (23) (2004),
4953-4966.

[6] K. Yakut, N. Alemdaroglu, I. Kotciolglu, C. Celik, Expérimental investigation of

thermal résistance of a heat sink with hexagonal fins, Applied Thermal Engineering
26 (17-18)(2006), 2262-2271.
[7] Y. Chen, M. Groll, R. Martz, Y. F. Maydanik, S.V. Vershinin, Steady-state and
transient performance of a miniature loop heat pipe, International Journal of ther­
mal Sciences 45 (11)(2006), 1084-1090.
[8] S. Hong, D. R. Herling, Open-cell aluminium foams filled with phase change ma-
terials as compact heat sinks, Scripta Materialia 55 (10)(2006), 887-890.

96
 97

[9] R. Chein, J. Chuang, Expérimental microchannel heat sink performance studies

using nanofluids, International Journal of Thermal Sciences 46 (1)(2007), 57-66.
[10] K. Boomsma, D. Poulikakos, F. Zwick, Métal foams as compact high performance
heat exchangers, Mechanics of Materials 35 (12)(2003), 1161-1176.
[11] G. Hetsroni, M. Gurevich, R. Rozenblit, Sintered porous médium heat sink for
cooling of high-power mini-devices, International Journal of heat and fiuid Flow
27 (2)(2006), 259-266.
[12] D. P. Haack, K.R. Butcher, T. Kim, T.J. Lu, Novel lightweight métal foam ex­
changers,American Society of Mechanical Engineers, Process Industries Division
(Publication) PID 6 (2001), 141 -147.
[13] A. Bejan, Designed porous média : maximal heat transfer density at decreasing
length scales, International Journal of Heat and Mass Transfer 47 (14-16)(2004),
3073-3083.
[14] Y.S. Muzychka, Constructal multi-scale design of compact micro-tube heat
sinks ans heat exachangers, International Journal of thermal Sciences (2006),
doi :10.1016/j.ijthermalsci.2006.05.002, article sous presse.
[15] A. Bejan, Shape and Structure, from Engineering to Nature, Cambridge University
Press, 2000.
[16] X. Wei and Y. Joshi, Optimization Study of Stacked Micro-Channel Sinks for
Micro-Electronic Cooling, IEEE Transactions on Compoents and Packaging Tech­
nologies 26 (1)(2003), 55-61.
[17] K. Jeevan, G.A. Quadir, K.N. Seetharamu, LA. Azid, Z.A. Zainal, Optimisation of
thermal résistance of stacked micro-channel using genetic algorithms, International
Journal for Numerical Methods in Heat and Fluid Flow 15 (1)(2005), 27-42.
[18] K. Foli, T. Okabe, M. Olhofer, Y. Jin, B. Sendhoff, Optimisation of micro heat
exchanger : CFD, analytical approach and multi-objective evolutionary algorithms,
International Journal of Heat and Mass Transfer 49 (5-6) (2006), 1090-1099.
[19] P. Jiang, M. Fan, G. Si, Z. Ren, Thermal-hydraulic performance of small scale
micro-channel and porous-media heat-exchangers, International Journal or Heat
and Mass Transfer 44 (5)2001), 1039-1051.
[20] A. Bejan, Convection Heat Transfer, 2d édition, Wiley, 1999.
[21] C. Hicks, A Genetic Algorithm tool for optimising cellular or functional layouts
in the capital goods industry, International Journal of Production Economies 104
(2) (2006), 598-614.
 98

[22] K. Shintania, A. Imaib, E. Nishimurab, S. Papadimitriouc, The container ship-

ping network design problem with empty container repositioning, Transportation
Research Part E : Logistics and Transportation Review 43 (1)(2007), 39-59.

[23] T. Fuhner, A. Erdmann, Improved mask and source représentations for automatic
optimization of lithographie process conditions using a genetic algorithm, Procee-
dings of SPIE - The International Society for Optical Engineering 5754 (1)(2005),
415-426.

[24] R. Le Riche, R.T. Haftka, Optimization of stacking séquence design for buckling
load maximisation by genetic algorithms, AIAA Journal 31 (5) (1993), 951-956.
[25] Z. Giirdal, R.T. Haftka, S. Nagendra, Genetic algorithm for the design of laminated
composite panels, SAMPE 30 (3)(1994), 29-35.
[26] G. Soremekun, Z. Gurdal, R.T. Haftka, L.T. Watson, Composite laminate design
optimisation by genetic algorithm with gêneralized elitist sélection, Computers and
Structures 79 (2)(2001), 131-143.
[27] R. Le Riche, R.T. Haftka, Improved genetic algorithm for minimum thickness
composite laminate design, Composite Engineering 5 (2)(1995), 143-161.
[28] V. G. Gantovnik, Z. Giirdal, L.T. Watson, A genetic algorithm with memory for
optimal design of laminated sandwich composite panels, Composite Structures 58
(4) (2002), 513-520.
[29] L. Valdevita, A. Pantanob, H. A. Stonea, A.G. Evansc, Optimal active cooling per­
formance of met allie sandwich panels with prismatic cores, International Journal
of Heat and Mass Transfer 49, (21-22) (2006), 3819-3830.
[30] T. Dias Jr., L. F. Milanezb, Optimal location of heat sources on a vertical wall
with natural convection through genetic algorithms, International Journal of Heat
and Mass Transfer 49 (13-14)(2006), 2090-2096.
[31] R. Selbasa, O. Kizilkana, M. Reppichb, Chemical Engineering and Processing 45
(4) (2006), 268-275.
[32] F. Girard, Optimisation de laminés en utilisant un algorithme génétique, Mémoire
de maîtrise, Faculté des sciences et de génie, Université Laval, 2006.
[33] L. Wanga, L. Zhangb, D. Zhenga, An effective hybrid genetic algorithm for flow
shop scheduling with limited buffers, Computers and Opérations Research 33
(10)(2006), 2960-2971.
[34] D. Neild, A. Bejan, Convection in Porous Media, Springer, New-York, 1999.
 99

[35] E. 0 . Holdzbecher, Modeling Density-Driven Flow in Porous Media, Springer, Ber­

lin, 1998.
[36] H. Darcy, Les Fontaines Publiques de la Ville de Dijon Victor, Dalmont, Paris,
1856.
[37] F. P. Incropera, D. P. DeWitt, Fundamentals of Heat and Mass Transfert, 5th Ed.,
John Wiley and Sons, 2002.
[38] D. C. Pham, S. Torquato, Strong-contrast expansions for the effective conductivity
of isotropic multiphase composites, Journal of Applied Physics 94 (10) (2003), 6591-
6602.
[39] B. R. Munson, D. F. Young, Théodore H. Okiishi, Fundamental of Fluid Mechanics,
4th Ed., John Wiley & Sons inc, New-York, 2002.
[40] www.metalprices.com [08/2006]
[41] J. H. Ferziger, Milovan Peric, Computational Methods for fluid Dynamics, 3rd Ed.,
Springer, 2002.
[42] Fluent 6.1.22, Fluent user's Guide, Fluent Incorporated,www.fluent.com, 1998.
[43] Gambit 2.2.30, www.fluent.com.
[44] A. Chipperfield, P. Fleming, H. Pohlheim, C. Fonseca, Genetic Algorithm Tool-
box For Use with Matlab, IEE Colloquium on Applied Control Techniques Using
Matlab 14 (1995),10/-10-4.
[45] J.J. Grefenstette, Optimization of control parameters for genetic algorithms, IEEE
Transactions on Systems, Man, and Cybernetics 16 (1)(1986), 122-128.
[46] J. Chen, X. Yang, Optimal parameter estimation for Muskingum model based
on Gray-encoded genetic algorithm, Communications in Nonlinear Science and
Numerical Simulation (2005),doi :10.1016/j.cnsns.2005.06.005, article sous presse.
[47] D.Whitley, The Genitor Algorithm and Sélection Pressure : Why Rank-Based Al­
location of Reproductive Trials is Best, Third International Conférence on Genetic
Algorithms, San Mateo, California, USA, Morgan Kaufmann Publishers 1989, 116-
121.
[48] J. E. Baker, Reducing Bias and Inefnciency in the Sélection Algorithm, Proceedings
of the Second International Conférence on Genetic Algorithms and their Applica­
tion, Hillsdale, New Jersey, USA, Lawrence Erlbaum Associates 1987, 14-21.
[49] F. Bobaru, S. Rachakonda, Optimal shape profiles for cooling fins of high and
low conductivity, International Journal of Heat ans Mass Transfer 47 (23) (2004),
4953-4966.
[50] www.zalmanusa.com [09/2006].

[51] A. A. Mohammadein, M. F. El-Amin, Thermal radiation effects on power-law fluids

over a horizontal plate embedded in a porous médium, International Communica­
tions in Heat and Mass Transfer 27 (7)(2000), 1025-1035.

[52] S. M. Al-Harbi, Numerical study of natural convection heat transfer with variable
viscosity and thermal radiation from a cône and wedge in porous média, Applied
Mathematics and Computation 170 (1)(2005), 64-75.

[53] I. A. Badruddina, Z.A. Zainala, P.A. A.Narayanab, K.N. Seetharamub, Numerical

analysis of convection conduction and radiation using a non-equilibrium model in
a square porous cavity, International Journal of Thermal Sciences 46 (1)(2007),
20-29.
[54] A. C. Baytas, I. Pop, Free convection in a square porous cavity using a thermal
nonequilibrium model, International Journal of Thermal Sciences 41 (9) (2002),
861-870.

[55] N. H. Saeida, A. A. Mohamadb, Periodic free convection from a vertical plate in

a saturated porous médium, non-equilibrium model, International Journal of Heat
and Mass Transfer 48 (18) (2005), 3855-3863.
[56] M. F. El-Amin, Double dispersion effects on natural convection heat and mass
transfer in non-Darcy porous médium, Applied Mathematics and Computation
156 (1)(2004), 1-17.

[57] J.K. Cochran, S.M. Horng, J.W. Fowler, A mufti-population genetic algorithm to
solve multi-objective scheduling problems for parallel machines, Computers and
opérations research 30 (7) (2003), 1087 -1102.
[58] N.S. Mera, L. Elliott, D.B. Ingham, A mufti-population genetic algorithm approach
for solving ill-posed problems, Computational mechanics 33 (4) (2004), 254-262
 Annexe A

Exemple de problème dimensionnel

À titre d'exemple dimensionnel étudions le stratifié présenté à la figure 5.6. Cette

échangeur opère à RaD2 = 10 9 , et possède les caractéristiques suivantes : L = 1.0,
N = 10, D = 0.001, Ra = 10 15 , fmax = 3.96E - 4. Considérons l'air comme fluide
refroidisssant. Dans le cas où on limite la différence de température maximale T — T^ à
200 degrés, on obtient un échangeur ayant une hauteur de H = 3.4 m pouvant supporter
un flux atteignant q" = 5.1 kW/m2. On constate qu'en raison de la taille du refroidis-
seur, celui-ci ne peut être utilisé dans le refroidissement de composantes électroniques.
Cette configuration se prête mieux au refroidissement de système de grande dimension,
comme ceux que l'on retrouvent par exemple dans des centrales thermiques. Il est à
noter que lors des simulations présentées dans ce document, la dimension des pores est
considérée fixe (D = 0.001). D'après les choix de quantités de référence (éq. 2.19-2.23)
employées lors de l'adimensionalisation de l'équation d'énergie éq. 2.25 et l'équation
de Darcy éq. 2.24, le système thermofluide est en fait régie par la multiplications des

101
paramètres RaD2. Il est donc possible de considérer des résultats d'optimisation cor­
respondant à des configurations poreuses ayant des valeurs différentes de D auquelles
on associe un nombre de Rayleigh approprié, en autant que la combinaison RaD2 soit
la même. Ainsi, dans le contexte du refroidissement de composantes électroniques, où
la taille des échangeurs est de l'ordre de 10 à 50 cm, il faudrait effectuer des simulations
pour des nombre de Rayleigh de l'ordre 1012 — 1013 lorsque (D = 0.001).
 Annexe B

Script de l'algorithme génétique

Cette annexe comprend une copie du script employé pour résoudre le problème
d'optimisation par algorithme génétique. On retrouve d'abord un résumé des principales
fonctions (section B.l). Ensuite, un exemple de la séquence des opérateurs utilisés dans
l'AG est présenté (section B.2). À la section B.3 chacun les opérateurs sont détaillés
selon leurs apparitions dans la routine d'optimisation.

B.l R é s u m é du contenu de l'AG

'A Projet de Maîtrise:

103
/, Optimisation d'un stratifié poreux
*/, soumis à un flux thermique convection naturelle
•/. Date: 15.12.2006
'/. Université Laval, Québec
y.
'/, Ensemble de fonction constituant L'algorithme génétique.
%
y.
y.
'/, Création de population initiale
%
'/. crtbp - Crée une population initiale en binaire
'/, crtbase - Crée un vecteur de la base
y.
°/„ Décoder la population
t
'/, bs2rv.m - Transforme une chaîne de caractère
"L en vecteur de valeurs réelles
'/. rv2bs.m - Traduit la phénotype en chromosome
'/, enigma.m - Extrait la représentation binaire en code gray
'/. d'un nombre réel
y.
% Corriger les individus
t
'/, correctgene.m - Cette fonction modifie les gènes d'un chromosome
'/, Dépourvus de sens physique.
%
'/. troque.m - Cette fonction tronque les phénotypes qui
'/. possèdent une couche vide au sein de leur
% architecture.
y.
7t Fonction objectif
%
'/, naturelle.m - Résout le système thermofluide formé
y, formé d'un échangeur poreux en convection
*/■ naturelle
'/, pénalité.m - Modifie la valeur de fonction objectif si les
Z contraintes ne sont pas respectées
% matériaux.m - Banque de matériaux disponible
'/. lors de l'optimisation
'/, propriété.m - Charge les propriétés du matériau en lien
'/. avec la banque de matériaux
% conductivite.m - Calcule la conductivité équivalente
% perméabilité.m - Calucle la perméabilité de la structure
y.
% Classement des individus
%
'/. ranking.m - Assigne un indice selon la performance
y.
y. Sélection and reinsertion
y,
'/, reins.m - uniform random and fitness-based reinsertion
'/, sélect.m - high-level sélection routine
'/, sus.m - stochastic uni versai sampling
X Opérateurs de Mutation
y,
"/. mut.m - Effectue une mutation discrète
y.
°/. Opérateur de croisement
%
y, xovmp.m - croisement multipoint
y.
X Recherche locale
t
*/, locale.m - Effectue une recherche locale
y.
'/. Banque de données
%
Z histoire.m - Effectue à la gestion des phénotypes évalués X
% Fonction utile
%
'/, trapèze.m - Effectue un intégrale d'une fonction grâce à
'/. méthode des trapèzes
°/, tridiagC - Effectue la résolution d'un système matriciel
% tridiagonal
'/. Fin

B.2 Boucle principale

'/. genopt.m - optimisation avec contraintes

y.
y. Cette fonction effectue l'optimisation des couches d'un
'/. stratifier poreux afin de minimiser la température maximale
'/. atteinte à la paroi où un flux de chaleur est impose.
'/. Les variables de l'optisation sont la porosité et
'/, le type de matériau qui compose chaque couche.
y.
'/, Auteur : Charles Villemure
'/, Historique: 18.07.05 : création de la fonction

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxy.y.xxxxxx

clear
clc
I; i c
graphique=0;

X initialisations des paramètres

[strat.algo] = initstrucQ; X structure qui résume les résultats
'/, strat — > propriétés du stratifié opt.
X algo — > propriétés de l'AG
'/, boucles sur différentes valeurs du paramètre alpha
X de la fonction objectif modifiée, opt. masse et coût sous cont. de temp.
 106

for ialpha=2:2
'/, dimension des populations
NIND = 3; '/, Nombre d'individus par sub-population
algo.nind=NIND;
MAXGEN = 5; y Nombre maximal de générations
algo.maxgen=MAXGEN;
MAXarret = 250; y Nombre génération sans changement
algo.maxarret=MAXarret;

% Paramètres du Stratifié
NSTRAT=10; '/. Nombre de couche stratifié
algo.nstrat=NSTRAT;
strat.nstrat=NSTRAT;

■/. Plage des différentes longueurs du stratifié ( boucle )

longueur=[0.25];

7, Boucle sur le nombre de Rayleigh

% Plage des différents nombres de Rayleigh a couvrir ( boucle )
Rayleigh=[lel5];

'/. Les matériaux

7. Banque de Matériaux
nbMat=4; "/„ nombre de matériaux a considérer : 4 ou 8
algo.nbmat=nbMat;

'/, Les contraintes de masse ( boucle)

'/, fraction de la masse sans contrainte massex*masseScont = masseMax
massex=[0.5];

'/, Les contraintes de coût ( boucle)

"/, fraction du coût sans contrainte coutx*coutScont = coutMax
coutx= [] ;

X détermine si une double contrainte est imposée Masse et coût

'/. 1 : active le combo '/, Scont + %cont
'/, 2 : active le combo '/, cont + '/, cont même degré
Z 3 : active le combo % cont + 7% cont couplé
'/, autre : non actif
combo=l; '/. détermine quel sorte de test sera effectué ( 1, 2 ou 3 )
strat.combo=combo;

°/. Les contraintes de température

% fraction de la température sans contrainte tempx*TempScont - TMaxlimite
'/. pour que ce type de simulations soit possible, il faut que les <
'/, vecteur massex et coutx soient vides.
tempx= [] ;
'/, identifier le nombre de simulations à réaliser
nbDecontrainte=length(massex)+length(coutx)+length(tempx);

if combo==3 y les vecteurs massex et coutx sont couplés

nbDecontrainte=nbDecontrainte-length(coutx);
end
if nbDecontrainte == 0 '/. dans le cas d'une optimisation non contrainte
nbDecontrainte=l;
end
'/, boucle globale sur les contraintes
for icont = l:nbDecontrainte;

% Le nom de la simulation
nom_simulation=sprintf( 'maptoutecontrainte05_y,d' ,icont) ;
X gestion des répertoires ( création d'un nouveau dossier )
nom_simulation_rep=repertoire(nom_simulation);

'/,% boucle sur les différentes longueurs d'échangeur '/.

for ilong=l:length(longueur)
'/, Longueur de l'échangeur
L=longueur(ilong);
strat.1=L;

'/. Porosité maximale : correspondant a une couche vide

poroslté_max=l.0;
algo.porositemax=porosite_max;

'/, Gestion des bornes des matériaux

MatINF=l; '/. matériaux disponibles : Borne inférieure
MatSUP=MatINF+ (nbMat-1); '/. matériaux disponibles : Borne supérieure

if nbMat==4
PRECI_matID = 2; '/. Précision de la représentation binaire
'/. pour la variable : matlD
elseif nbMat==8
PRECI_matID = 3; '/, Précision de la représentation binaire
'/, pour la variable : matlD
end
algo.precimatid=PRECI_matID;

"/, Possibilité de plusieurs matériaux

Borne_MatINF=MatINF;
Borne_MatSUP=MatSUP;

'/, boucle sur le nombre de Rayleigh

for iRay=l:length(Rayleigh)
Ra=Rayleigh(iRay);
strat.ra=Ra;

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
TIXIXIXIXIXIXIXIX Algorithme Génétique XXXXXXXXXXXXXXXXXXXXX
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
■/, Les contraintes
X Optimiser Tmax sous des contraintes de [masseMax, coutMax]

"/, déterminer les valeurs de la masse et du coût maximal selon Ra et L

X récupère les valeurs des propriétés de l'échangeur pour un Ra et un L
X à partir des résultats sans contraintes
[TmaxScont.masseScont,coutScont]=PropSansCont(Ra,L);
strat.tmaxscont=TmaxScont;
strat.massescont=masseScont;
strat.coutscont=coutScont;
7. pourcentage de la contrainte par rapport au résultat sans contraintes
rm=0; 7. s'adresse à la masse
rc=0; '/. s'adresse au coût
rt=0; 7. s'adressse à la température

7. déterminer les contraintes à imposer

7. contrainte de masse
if icont <= length(massex) kk length(massex)~=0
masseMax=massex(icont)*masseScont;
rm=massex(icont);
else
masseMax=0;
rm=0;
end

7. contrainte de coût
if icont > length(massex) kk length(coutx)~=0
coutMax=coutx(icont-length(massex))*coutScont;
rc=coutx(icont-length(massex));
else
coutMax=0;
rc=0;
end

7. contrainte sur la température

if icont > (length(massex)+length(coutx)) kk length(tempx)~=0
TMaxlimite=tempx(icont-length(massex)-length(coutx))*TmaxScont;
rt=tempx(icont-length(massex)-length(coutx));
alpha=0.01*(i alpha-1);
else
TMaxlimite=0;
rt=0;
alpha=0;
end

7. contrainte de masse et de coût combiné : approche 1

7. pourcentage de SCO identique pour une des deux contraintes (masse et coût)
7. l'autre contrainte correspond à 507. de la propriété SCO
if combo==l
7. récupère les valeurs des propriétés de l'échangeur pour un Ra et un L
7. à partir des résultats avec contraintes
[TmaxAveccont.masseAveccont,coutAveccont]=PropAvecCont(Ra,L,rm,rc);
if rm"=0.0
coutMax=0.5*coutAveccont;
elseif rc~=0.0
masseMax=0.5*masseAveccont;
end
end 7. if combo==l

7. contrainte de masse et de coût combiné : approche 2

7. même pourcentage de SCO masse que SCO coût
if combo==2
 '/. récupère les valeurs des propriétés de l'échangeur pour un Ra et un L
'/, à partir des résultats avec contraintes
if rm~=0.0
coutMax=rm*coutScont;
elseif rc~=0.0
masseMax=rc*masseScont;
ond
end '/, if combo==l

'/, contrainte de masse et de coût combiné : approche 3

'/. deux pourcentages différents de SCO sur la masse et le coût ( couplés )
if combo==3
'/, récupère les valeurs des propriétés de l'échangeur pour un Ra et un L
'/. à partir des résultats avec contraintes
masseMax=massex(icont)*masseScont;
coutMax=coutx(icont)*coutScont;
rm=massex(icont);
rc=coutx(icont);
end % if combo==3

'/, mise à jour

strat.massemax = masseMax;
strat.coutmax = coutMax;
strat.tmaxlimite = TMaxlimite;
strat.rm = rm;
strat.rc = rc;
strat.rt = rt;
strat.alpha = alpha;

y, loger les limites dans un vecteur masse, coût, temps

contraints=[masseMax,coutMax,TMaxlimite,alpha];

Z paramètre d'arrêt de la boucle : lorsque arrêt == 1 — > stop

arret=0;

'/, initialiser la première centaine ( paramètre du graphique )

uncentaine = 0;

'/. identifier s'il y a un problème dans les rouages de l'AG

problème = 0;

°/, indicateur du nombre de fois que la valeur de ObjV demeure inchangées

y, lorsque inchangée atteint le critère d'arrêt — > stop
inchangé = 1 ;

'/, indicateur du nombre de fois que l'amélioration locale est bénéfiques

stumieux = 0;
algo.stumieux=stumieux;

% Paramètres du stratifié
NSTRAT=10; '/. Nombre de couche stratifié

"/, connaître les propriétés des matériaux disponibles

[tableDdensite.tableDk.tableDmat.tableDcout,nbDmat]= matériaux ( 1 ) ;
7. Paramètres de l'algorithme génétique
GGAP = 0.8; 7. Génération gap: Le nombre d'individus crée
algo.ggap=GGAP;
Pmut = 0.05; 7. Probabilité de mutation
al go.pmut=Pmut;
Ncrois = 5 ; 7. nombre de coupe dans le croisement
algo.ncrois=Ncrois;
Pcrois = 0.8 ; 7. Prababilité de croisement
algo.pcrois=Pcrois;
Genlocale = 1; 7. nombre de générations entre deux opt. local
algo.genlocale=Genlocale;

7. derniers matériaux dans la liste

MatDernier=4; 7. derniers matériaux possible dans la liste
% après quoi on obtient un matériau vide
algo.matdernier=MatDernier;

7. Spécification des variables

NVAR_matID = NSTRAT; 7. Nombre de variables
NVAR_poro = NSTRAT; 7. Nombre de variables

7.PRECI_matID = 2; 7. Précision de la représentation

7 binaire pour la variable : matlD
DOM_matID=[Borne_MatINF;Borne_MatSUP] ;7. Domaine de la variable matlD
PRECI_poro = 5; 7. Précision de la représentation
algo.preciporo=PRECI_poro; 7. binaire pour la variable : poro
D0M_poro = [0;l]j 7. Domaine de la variable poro

7. Spécifier les routines employées

SEL_F = 'sus'; 7. Type de sélection
X0V_F = 'recdis'; 7. Recombinaison en fonction des individus
MUT_F = 'mutbga' ; 7. Type de Mutation
0BJ_F = 'naturelle'; 7. nom de la fonction objective

7. Initialisation du Compteur et de la liste des meilleurs

Meilleur = NaN*ones(MAXGEN,l) ; 7. meilleurs individus
algo.meilleur=Meilleur; 7. de chaque génération
Moyenne = NaN*ones(MAXGEN, 1) ; 7. Moyenne des performances
algo.moyenne=Moyenne; 7. des inds de chaque gen.
gen = 0; 7. compteur de génération

7. base de donnée initiale (nulle)

[ObjVhist,Tmaxhist,Dbase]= historique(1,0,-1,NSTRAT,Ra,L.masseMax,coutMax);

7. début de A . G.

7. Construction du "champ descripteur"

FieldD = [rep( [PRECI_matID;DDM_matID;1;0;1;1],[l,NVAR_matID]),,
[PRECI_poro;DDM_poro ; 1 ; 0; 0; 0 ] , . . .
rep([PRECI_poro;D0M_poro;l;0;l;l],[l,NVAR_matID-l] )] ;

7. Traduction du binaire au réel

[traduction,Phen_possible] ■ enigma(FieldD);

7. Création de la population initiale

 111

Chrom = crtbp(NIND, (NVARjiatID*PRECI_matID + NVAR_poro*PRECI_poro));

7. Phenotype de la population initiale

Phen = bs2rv(Chrom,FieldD);

'/, Correction des gènes

[Phen,Chrom] = correctgene(Phen,Chrom,FieldD,...
PRECI_matID,PRECI_poro,NSTRAT,...
MatDernier,porosite_max);

'/, Tronquer le stratifié si l'on rencontre une couche a porosité maximale

[Chrom,Phen] = tronque(Chrom,Phen,PRECI_matID,PRECI_poro,...
porosite_max,NSTRAT,MatDernier,...
traduction,Phen_possible);

7. Évaluation de la population initiale

[DbjV.Dbase.T,v,Prop,strat,algo] = n a t u r e l l e ( s t r a t , a l g o , . . .
Phen,NSTRAT,Ra,Dbase,L,contraintes);
if graphique==l
7, Déterminer le meilleur individu et afficher la convergence
figure('Name',['Evolution des Meilleurs'],'NumberTitle',['off'])
title(['Le meilleur du la génération ']);
Moyenne(gen+l)=mean(ObjV);
[Meilleur(gen+1).indice] - min(ObjV);
allureMeilleur(gen+l,:)=Phen(indice,:) ;
hold on
plot((Meilleur),'ro');
plot( Moyenne,'b+');
hold off
xlabeK 'génération' ) ; ylabeK 'loglO(f (x)) ') ;
text(0.5,0.90,['Meilleur = ' , . . .
num2str(Meilleur(gen+l))],...
'Units','normalized');
text(0.5,0.85,['Génération = ' , num2str(gen)],'Units','normalized');
drawnou;
clf
end '/, if graphique==l
7. Boucle sur les Générations
while arrêt < 1 && probleme==0

7. Assigner une valeur de "fitnesse" a l'ensemble de la population

FitnV = ranking(ObjV);

7. Sélectionner les individus pour la reproduction

SelCh ■ sélect('sus', Chrom, FitnV, GGAP);

7. Amélioration locale à toute les 25 générations

if rem((gen+l)/Genlocale,2)==0

7. Correction des gènes

PhenRejet = bs2rv(SelCh,FieldD);
[PhenRejet,SelCh] = correctgene(PhenRejet,SelCh,FieldD,...
PRECI_matID,PRECI_poro,...
NSTRAT,MatDernier,porosite_max);
 '/, Tronquer le stratifié
'/, si l'on rencontre une couche à porosité maximale
[SelCh.PhenRejet] - tronque(SelCh,PhenRejet,PRECI_matID,...
PRECI_poro,porosite_max,...
NSTRAT.MatDernier);

PhenRejet = bs2rv(SelCh,FieldD);
[PhenRejet,SelCh,Dbase,stumieux,strat.algo] = locale(...
PhenRejet,SelCh,...
FieldD,...
PRECIjnatID,PRECI_poro,...
NSTRAT.MatDernier,porosite_max,...
Ra,Dbase,L,contraintes,stumieux,...
strat.algo);

end % amélioration locale

■/. Recombiner les individus (crossover)

SelCh = xovmpCSelCh, Pcrois, Ncrois, 1 ) ;

°/, opérer la mutation sur les rejetons

SelCh = mut(SelCh,Pmut);

"/. Correction des gènes

PhenRejet = bs2rv(SelCh,FieldD);
[PhenRejet,SelCh] = correctgene(PhenRejet,SelCh,FieldD,...
PRECI_matID,PRECI_poro,...
NSTRAT.MatDernier,porosite_max);

'/, Tronquer le stratifié si l'on rencontre une couche a porosité max.

[SelCh.PhenRejet] = tronque(SelCh,PhenRejet,PRECI_matID,...
PRECI_poro,porosite_max,...
NSTRAT.MatDernier);

'/, Évaluation des rejetons

[ObjVSel.Dbase.T,v.Prop,strat.algo] = naturelle(strat,algo,...
PhenRejet.NSTRAT.Ra,...
Dbase.L,contraintes);

'/. Réinsérer les rejetons dans la présente population

[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,DbjVSel);
Phenfinal = bs2rv(Chrom,FieldD);

y, incrémenter le compteur gen = gen+1;

display(gen);

y. Mise a jour de la visualisation et de notre " plus Meilleur "

if rem(gen/100,l)==0
uncentaine=uncentaine+l ;
end

if MAXGEN <100
topPlot=MAXGEN;
else
topPlot=100;
 end
[Meilleur(gen+l).indice] = min(ObjV);
Moyenne(gen+l)=mean(ObiV);
allureMellleur(gen+l,:)=Phenfinal(indice,:);

if (Meilleur(gen+l)>Meilleur(gen+D) ...
I isnan(Moyenne(gen+D) I isnan(Meilleur(gen+l))
probleme=l;
'/. Savegarde des donne en b i n a i r e dans nom_fichier
nom_fichier=sprintf (''/.sWERREUR.Ra'/M.Od.L^/.l. ld_'/,s.mat ' , . . .
nom_simulation_rep,Ra,L,datestr(now, 30));
save(nom_fichier);

clear Dbase Moyenne Meilleur T v Prop LeMeilleur

clear allureMeilleur
end % problème

if graphique==l
'/, graphique
hold on
plot( Meilleur(uncentaine*100+l:uncentaine*100+topPlot),'ro');
plot( Moyenne(uncentaine*100+l:uncentaine*100+topPlot),'b+');
hold off
x l a b e H ' g é n é r a t i o n ' ) ; y l a b e K ' (f (x)) ' ) ;
t i t l e ( [ ' L e m e i l l e u r de l a g é n é r a t i o n ']);
text(0.5,0.90,['Meilleur = ', ...
num2str(Meilleur(gen+l))],'Units','normalized');
t e x t ( 0 . 5 , 0 . 8 5 , ['Génération = ' , num2str(gen)],'Units','normalized');
drawnow;
c If-
end '/. if graphique==l

7, critère d'arrêt
if Meilleur(gen+1)== Meilleur(gen)
% incrémenter le compteur
inchanger=inchanger+l;
else
% remise du compteur à zéro
inchanger=l;
end
'/. arrêt si le critère de convergence est atteint if inchanger > MAXarret
arret=l;
end
'/„ arrêt si le nombre de générations maximales est atteint
if gen >= MAXGEN
arret=l;
end

7. affichage à l'écran
if rt"=0
gen Meilleur inchanger= [gen,Meilleur(gen+l)/1000,inchanger] ;
else
gen Meilleur inchanger=[gen,Meilleur(gen+l)*1000,inchanger] ;
end
display( gen Meilleur inchanger) ;
 114

end '/.boucle while

if probleme==0
% Déterminer l ' a l l u r e de l a d e r n i è r e génération
disp(indice);
dispOCa c e s t n o t r e "plus m e i l l e u r 1 " ) ;
display(indice);
d i s p C L a configuration optimale e s t : ' ) ;
for i=l:NVAR_matID
disp(tableDmat(Phenf i n a K i n d i c e , ( i ) ) ) ) ;
disp(Phenf i n a K i n d i c e , (NVAR_matID+i))) ;
end
dispCTempérature Maximale ' ) ;
disp(Meilleur(gen+l));
[LeMeilleur,Dbase,T,v,Prop,strat,algo] = naturelle(strat,algo,...
Phenf i n a K i n d i c e , : ) , . . .
NSTRAT.Ra.Dbase,. . .
L,contraintes,1) ;
algo.alluremeilleur=allureMeilleur;
algo.meilleur=Meilleur;
algo.moyenne=Moyenne;

strat.distmat=Phenfinal(indice,1:NSTRAT);
s t r a t . d i s t p o r o = P h e n f i n a K i n d i c e , (NSTRAT+1) :2*NSTRAT) ;
s t r a t . t m a x = Prop(3);
strat.masse = Prop(l);
s t r a t . c o u t = Prop(2);

'/, Savegarde des donne en binaire dans nom.fichier

nom_fichier=sprintf ('y„s\\Ray„4.0d_L=y.l.ld_y.s.mat',. . . .
nom_simulation_rep,Ra,L,datestr(now, 30));
save(nom_fichier);

clear Dbase Moyenne Meilleur T v Prop LeMeilleur

clear allureMeilleur

'/, Effacer les variables de grandes tailles génèrent lors de la simulation

% clear Dbase allureMeilleur T T_ancien Meilleur
end '/, problème ==0
end '/, boucle sur les Ra
end '/. boucle sur les L
end '/, icont boucle sur les contrainte
end Z iNSTRAT
'/, End of GA
toc
t=toc;
•/. Fin
B.3 Fonctions de l'AG

B.3.1 Banque de matériaux

% matériaux.m - dresse liste des propriétés de différents matériaux

'/• Disponible lors de l'optimisation.
%
% Cette fonction dresse liste des propriétés de différents matériaux et les
'/, organisent dans un tableau.
•/.
Z Syntaxe:
% [tableDdensite,tableDk,tableDmat,tableDcout,nbDmat]= materiaux(CHX)
'/.
% Paramètre d'entrée:
%
'/. CHX Choix du type de valeurs retourné par la fonction
z si CHX = 1 : la fonction retourne la table des
1 densités, la table des matériaux, la table des
V. Conductivités ainsi que le nombre de matériaux.
1 Si CHX = 2 : la fonction retourne seulement
y, le nombre de types de matériaux.
y.
'/, Paramètre de sortie:
•/.
% tableDdensite - Vecteur contenant la valeur de la densité se
7. rattachant au matériau demandé
I,
7. par l e l ' i n d e n t i f i a n t matlD
%
'/. tableDk - Vecteur contenant la valeur de la
'/. conductivité thermique se rattachant au matériau
y. déterminé par le l'identifiant matlD
•/.
'/. tableDmat - Vecteur contenant les types de matériaux
y.
Z tableDmat - Vecteur contenant le coût de chaque matériel
Z
Z nbDmat - Le nombre de matériaux compris dans la table
y,
y.
Z Auteur : Charles Villemure
'/. 15.12.2006
Z
z************************+**+**********************************************
function [tableDdensite,tableDk,tableDmat.tableDcout,nbDmat]=materiaux(CHX)
y,**************************************************************************

Z Coûts normalisés par rapport au Fer: Al « 8.9, Cu = 18.7, Fe = 1,

'/. Brass = 11.6
Z
Z Propriété normalisée par rapport aux propriétés de l'air à 300 K.
i-i;
'/.l- ALUMINIUM
tableDmat(i)={'aluminium'};
tableDdensite(i)=2327.5;
tableDk(i)=90U.4;
tableDcout(i)=12.8;
i-i+ij

V.2- CUIVRE
tableDmat(i)={'cuivre'};
tableDdensite(i)=7691.6;
tableDk(i)=15247.2;
tableDcout(i)=18.9;
i-i+1;

7.3- ACIER
tableDmat(i)={'fer'};
tableDdensite(i)=6776.3 ;
tableDk(i)=3049.4;
tableDcout(i)=l;
i-i+i;
'/A- VIDE
tableDmat(i)={'vide'};
tableDdensite(i)=0;
tableDk(i)=0;
tableDcout(i)=0;
i-i+1I
if (length(tableDdensite)"=length(tableDk) &...
length(tableDmat)~=length(tableDk))
warning = ('trouble avec les indices');
end

nbDmat=length(tableDmat);

if CHX ==2
tableDmat='nulle';
tableDdensite='nulle';
tableDk='nulle';
tableDcout='nulle';
end

Z fin

B.3.2 Historique

'/, historique.m - gestion d'une base de donnée des individus testés

'/.
"/, Cette fonction vérifie si un individu possédant les mêmes
'/. caractéristiques a déjà été calculé auparavant.
'/, Si tel est le cas, la valeur de ObjV calculée antérieurement est fournie.
% Dans le cas contraire, la valeur de -1 est retournée.
 117

7. Syntaxe: [ObjVhist,Tmaxhist.Dbaseini] =
'/. historique (première, Phen, Dbase, NSTRAT, Ra, L .masseMax, coutMax)
7.
"/, Paramètre d'entrée:
y. •/.
'/. première - appeler la fonction pour la "première " fois afin de
7, Créer une basse de donnée initiale. Pour ce faire
'/, première = 1 , l'appel de la fonction
'/. Dbase = 0
y.
'/, Phen - table contenant les phénotypes des individus
*/. chaque ligne correspond a un individu
'/. chaque colonne correspond à une variable
'/, la première moitié des variables correspondent
'/, aux matériaux qui composent les couches
'/, la seconde moitié correspond
'/, à la porosité de chacune des couches
y.
'/, paramètre de sortie :
y.
'/, ObjVhist - valeur de référence sur laquelle repose l'optimisation du
'/, Design. Par exemple, celle-ci peut correspond
'/, a la valeur maximale de la
*/, température atteinte dans la distribution de température a
'/. laquelle deux pénalités sont ajoutées afin de
'/, répondre a des contraintes de coût et de poids du
'/. Stratifié. Si aucun individu ne correspond
'/, dans la base de donne, la valeur de -1 est retournée.
y.
% Tmaxhist - Valeur de la température maximale atteinte
'/ Dans le stratifié. Si aucun individu ne correspond
'/, dans la base de donne, la valeur de -1 est retournée.
y.
y.
'/. Dbaseini - base de données initiale comprenant une série de -1
y.
'/, Exemple d'une ligne de la base de données:
'/, Dbase= [matériaux des couches, porosité des couches,...
'/. Ra, L, mas seMax, coutMax, Obj V, Tmax]
'/, Auteur: Charles Villemure
•/. 20.12.2006

function [ObjVhist.Dbaseini] = h i s t o r i q u e ( p r e m i è r e , P h e n , n b v a r , . . .
Dbase, ObjV)

'/, initialisation de la base de données

if première == 1
'/, base de données initiale (nulle)
% Dbase= [matériaux des couches,porosité des couches, ...
7. Ra, L,masseMax,coutMax,ObjV,Tmax,]
clear Dbase
clear ObjVhist
Dbaseini(l,:) = [rep([-1] , [1,nbvar] ) ,-1] ;
 ObjVhist(l.l) = -1; '/. valeur de l'objectif
else
'/, Déterminer le nombre d'individus à vérifier
[m,n]=size(Phen) ;

'/. Création du tableau des individus traités

IndsPres=zeros(m,(nbvar)+l);
IndsPres(:,1 :(nbvar))=Phen;
IndsPresC:,nbvar+l)=ObjV;

'/. Déterminer les dimensions de la base de données

[tous,Ncaract]=size(Dbase);

'/, Vérifier si un individu de la base de données correspond

'/, à un individu traité

for idiv=l:m '/, boucle sur les individus présents

V, si aucun individu ne correspond dans la basse de données
ObjVhist(idiv.l)—1;
for j=l:tous '/, boucle sur la base de données
if( IndsPres(idiv,l:Ncaract-l)== Dbase(j,l:Ncaract-l) ) ;
ObjVhistCidiv,l)=Dbase(j,nbvar+l); '/, valeur de l'objectif
end 'h if == inds
end % boucle sur la base de donne
end '/, boucle sur les individus traites
end
'/. fin

B.3.3 Correspondance binaire-réel

'/. bs2rv.m - transforme une chaine de caractère

'/. en vecteur de valeurs réelles
y.
"/, Cette fonction décode les chromosomes de la forme binaire en une
'/, représentation réelle. Les chromosomes sont vus comme une chaîne
'/, de nombre binaire d'une certaine longueur encodée selon une
% structure numérale standard ou selon une strutrure Gray binaire.
'/. Ces chromosomes correspondent à des valeurs réelles comprises
°/0 dans un intervalle donné.
y
"/, Syntaxe: Phen = bs2rv(Chrom,FieldD)
y.
'/. Paramètres d'entée:
y.
'/, Chrom - Matrice contenant les chromosomes de la
y Population. Chaque linge correspond à la chaine
'/. de binaires associées à un individu.
'/, ( la première variable est à gauche
'/. et la dernière à droite )
y.
"/, FieldD - Matrice décrivant la longueur et comment décoder
'/. Chaque composante des chromosomes.
7,
'/• Selon la struture suivante:
'/.
7.
7. [len; (scalaire: nombre de bit de la variable)
7. lb; (scalaire: valeur de la borne inférieure du domaine)
7. ub; (scalaire: valeur de la borne supérieure du domaine)
7. code; (0=binaire I l=Gray)
7. scale; (0=arithmetique l l=logarithmique)
7. lbin; (0=inclule borne inf. I l=exclule borne inf. )
7. ubin] ; (0=lnclule borne sup. I l=exclule borne sup. )
7.
% où
7. len - vecteur rangé contenant le nombre de bit de chaque
7. des variables. La somme des longueurs
7. doit correspondre avec la longeur d'un individus.
7. lb - valeur de la borne inférieure du domaine
7. ub - valeur de la borne supérieure du domaine
7. code - Type de code : binaire ou Gray
7. scale - Type de discrétisation :
7. arithmétique ou logarithmique 7.
7. lbin - inclule ou exclure la borne inférieure.
7. ubin - inclule ou exclure la borne supérieure.
7.
'/.
% Output parameter:
7.
7. Phen - matrice contenant les valeurs
7. réelles des phénotypes
7.
7. Auteurs: Carlos Fonseca,Andrew Chipperfield
7. Adaptation: Charles Villemure
'/. 20.12.2006
7.
function Phen = bs2rv(Chrom,FieldD)
7.

7. nombre d'individus (Nind)

7. longueur des individus (Lind)
[Nind.Lind] = size(Chrom);

7. nombre de v a r i a b l e (Nvar)
[seven.Nvar] = s i z e ( F i e l d D ) ;

if seven "= 7
error('FieldD doit avoir 7 rangées.');
end

7. récupérer les propriétés de décodage

len = FieldD(l,:);
lb = FieldD(2,:);
ub = FieldDO, : ) ;
code = ~(TieldD(4, :));
scale - -(~FieldD(5,:));
lin = ~(~FieldD(6,:));
uin = ~(~FieldD(7,:));

% vérifié la concordance des propriétés,

if sum(len) "» Lind,
errorC'FieldD doit correspondre avec la longueur du chromosome ');
end

if ~all(lb(scale).*ub(scale)>0)
errorC'Échelle logarithmique ne peut contenir 0 dans le domaine');
end

7. décoder les chromosomes

Phen = zéros(Nind,Nvar);

lf = cumsum(len) ; '/, endroit ou couper les gènes

11 = cumsum([1 len]);
Prec = .5 ." len; 7. longueur des intervalles qui divise le domaine

logsgn = slgndb(scale)) ;
lb(scale) = log( abs(lb(scale)) ) ;
ub(scale) = log( abs(ub(scale)) ) ;
delta = ub - lb; '/, intervalle

Prec = .5 .* len;
num = ("lin) .* Prec;
den ■ (lin + uin - 1) .* Prec;

for i = l:Nvar,
idx - li(i) :lf (i) ;
if code(i) '/. code Gray
Chrom(:,idx)=rem(cumsum(Chrom(:,idx)')',2);
end
Phen(:,i) = Chrom(:,idx) * [ (.5)."(1:len(i))' ] ;
Phen(:,i) = lb(i) + delta(i) * (Phen(:,i) + num(i)) ./ (1 - den(i));
end

expand = ones(Nind,1);
if any(scale)
Phen(;,scale) = logsgn(expand,:) .* exp(Phen(:,scale));
end
7. FIN

B.3.4 Correspondance réel-binaire

7. rv2bs.m - cette fonction traduit la phénotype en chromosome

'/. à l'aide des propriétés binaires 7,
7. du gène, contenu dans le tableau FieldD.
7.
7. Syntaxe: [ Chrom ] = rv2bs( Phen, FieldD, ...
7. PRECI_matID, PRECI_poro, ...
7. NSTRAT, MatDernier)
 Paramètre d'entrée:

Phen table contenant les phénotypes des individus

chaque ligne correspond a un individu
chaque colonne correspond a une variable
la premier moitié des variables correspond au
matériaux qui composent les couches
la seconde moitié correspond a la porosité
de chacune des couches

FieldD caractéristique de la représentation

binaire d'un individu

y.
•/ ( une colonne pour chaque individu )
i
7, PRECI_matID Précision en nombre de bit
'/. pour le type de matériel 7.
7. PRECI_poro Précision en nombre de bit
7. pour la porosité 7.
7. NSTRAT Nombre de couches dans le stratifié
'/.
7. MatDernier indice du dernier matériau possible
7,
7. Paramètre de sortie:
y,
7. Chrom Matrice contenant la description des gènes
7. sous la forme de chromosomes chaque ligne
7. Correspondant a un individu.

7. Auteur : Charles Villemure

7. Historique: 09.07.06 : création de la fonction.
7.
7t******************************* *******************************************
function [ Chrom ] = rv2bs( Phen, FieldD, ...
PRECI_matID, PRECI_poro, ...
NSTRAT, MatDernier)
7t**************************************************************************
7. construire la table de traduction
[traduction,Phen_possible] = enigma(FieldD);

7. nombre de porosité possible

NbDePoro=length(Phen_possible(:,2));

7. dimension du tableau des phénotypes

[m,n]=size(Phen);

7. Construction du chromosome
Chrom=zeros(m,NSTRAT*(PRECI_matID+PRECI_poro));

for inds=l:m 7. boucle sur les individus

boutdechromo=zeros(l,PRECI_matID);
for igene=l : NSTRAT 7. boucle sur les gènes de matériaux
indice=find(Phen_possible(:,l)==Phen(inds,igene));
 if isempty(indice)
display(indice);
display(Phen(inds,igene) ) ;
pause(15)
save('Erreur');
end
boutdechromoC:,:)=traduction(1,indice,l:PRECI_matID);
ChromCinds, (PRECI_matID*(igene-l)+l) :PRECI_matID*(igene)) = . . .
boutdechromo;

end "/, for igene=l:n boucle sur les gènes

'A premier gène de porosité

boutdechrorao=zeros(1,PRECI_poro);
igene=NSTRAT+l;
indice=find(Phen_possibleC:,3)==Phen(inds,igene));
boutdechromoC:,:)=traduction(3,indice,1:PRECI_poro);
ChromCinds,CNSTRAT*PRECI_matID...
+ PRECI_poro*(igene-NSTRAT-l)+l) : ...
NSTRAT*PRECI_matID+PRECI_poro*(igene-NSTRAT))...
= boutdechromo;

for igene=(NSTRAT+2) :2*NSTRAT '/, boucle sur les gènes de porosité

indice=find(Phen_possible(:,2)==Phen(inds,igene));
boutdechromoC:,:)=traductionC2,indice,1:PRECI_poro);
ChromCinds,(NSTRAT*PRECI_matID...
+ PRECI_poro*Cigene-NSTRAT-l)+l) : ...
NSTRAT*PRECI_matID+PRECI_poro*Cigene-NSTRAT))...
= boutdechromo;
end '/. for igene=l:n boucle sur les gènes
end V, for inds=l:m boucle sur les individus

tester=0;
if tester==l
Phen_test=bs2rvCChrom,FieldD);
pareil=Phen_test-Phen;

for i=l;m
for j"i:n
if Phen_test~=0
displayCPhen_testCi,j));
end
end
end
end
'/.fin

B.3.5 Création de la population initiale

% crtbp.m - crée une population initiale en binaire

y.
% Cette fonction crée une population initiale en binaire
'/, selon une taille donnée et une structure particulière.
Z
'/. Syntaxe: [Chrom Lind BaseV] = crtbp(Nind, Lind, Base)
'/.
'/. Paramètres d'entée:
•/.
% Nind - Le nombre d'individus ou un vecteur rangé contenant le
'/. Nombre d'individus et leur longueur.
% Lind - Longueur du chromosome qui représente les individus
7, Base - Un scalaire ou un vecteur contenant la base des foyers
'/, Des gènes.
•/.
'/. Output Parameters:
•/.
'/. Chrom - Une matrice contenant des valeurs aléatoires
'/. de chromosome. Chaque individu est représenté
'/, sur une ligne. Matrice de dimension [Nind.Lind].
% Lind - Longueur du chromosome qui représente les individus
'/, BaseVec - Un vecteur qui correspond à la base des foyers des gènes
'/. associés à un chromosome.
*/.
'/, Auteur: Andrew Chipperfield
'/, Adaptation: Charles Villemure
"/. 20.12.2005
y,**************************************************************************
function [Chrom, Lind, BaseV] = crtbp(Nind, Lind, Base)
y^********************************************* ************ *****************
nargs = nargin ;

if nargs >« 1, [mN, nN] = size(Nind) ; end

if nargs >= 2, [mL, nL] = size(Lind) ; end
if nargs == 3, [mB, nB] = size(Base) ; end

if nN — 2
if (nargs = = 1 )
Lind = Nind(2) ; Nind = Nind(l) ; BaseV = crtbase(Lind) ;
elseif (nargs == 2 & nL == 1)
BaseV = crtbase(Nind(2),Lind) ; Lind = Nind(2) ; Nind = Nind(l) ;
elseif (nargs == 2 & nL > 1)
if Lind ~= length(Lind), errorCLind et la Base ne concordent pas'); end
BaseV = Lind ; Lind = Nind(2) ; Nind = Nind(l) ;
end
elseif nN == 1
if nargs == 2
if nL == 1, BaseV = crtbase(Lind) ;
else, BaseV = Lind ; Lind = nL ; end
elseif nargs == 3
if nB == 1, BaseV = crtbase(Lind,Base) ;
elseif nB ~= Lind, errorCLind et la Base ne concordent pas') ;
else BaseV = Base ; end
end
else
error('Erreur sur les paramètres en entrée') ;
end
"/, Crée une structure de chromosomes aléatoires de longueur [Nind.Lind]
/, La base des chromosomes est donnée par le vecteur ligne basse des foyers
Chrom = floor(rand(Nind,Lind).*BaseV(ones(Nind,1),:)) ;

7. fin

B.3.6 Correction des individus

correctgene.m - corrige la séquence de gènes

Cette fonction modifie les gênes d'un chromosome

afin que ceux-ci puissent demeurer dans la plage
souhaitée et conservent un sens physique.

Syntaxe: [PhenMOD] = correctgeneC Phen.Chrom.FieldD,...

PRECI_matID,PRECI_poro,
NSTRAT.MatDernier)

Paramètre d'entrée:3

Phen Table contenant les phénotypes des individus

chaque ligne correspond à un individu
chaque colonne correspond a une variable
la première moitié des variables correspond aux
matériaux qui composent les couches
la seconde moitié correspond a la porosité
de chacune des couches

Chrom Matrice contenant la description des gènes

sous la forme de chromosomes chaque ligne
Correspondant à un individu.

FieldD aracateristique de la représentation

binaire d'un individu
( une colonne pour chaque individu )

PRECImatlD Précision en nombre de bit

pour le type de matériau

PRECI_poro Précision en nombre de bit

pour la porosité

NSTRAT Nombre de couches dans le stratifié

MatDernier indice du dernier matériau possible

Paramètre de sortie:

Phen Phénotypes modifiés

Chrom Chromosome modifié

Auteur Charles Villemure

X 20.12.2006
%**************************************************************************
function [Phen,Chrom] = correctgene(Phen,Chrom,FieldD,...
PRECI_matID,PRECI_poro,...
NSTRAT,MatDernier,porosite_max)
X**************************************************************************
'/. construire la table de traduction
[traduction,Phen_possible] = enigma(FieldD);

X nombre de porosité possible

NbDePoro=length(Phen_possible(:,2));

X dimension du tableau des phénotypes

[m,n]=size(Phen);

'/, dimension du tableau des chromosomes

[mc,nc]=size(Chrom);

X copie de l'originale a modifier

PhenM0D=Phen;

for i = l:m X boucle sur les individus

'/, deuxième couche a une porosité maximale
paschanger=0;
if NSTRAT_=1
while paschanger==0
if (Phen(i,NSTRAT+2) == porosite_max )
X nouveau phenotype alleatoire
newPoroI =ceil(rand(l)*(NbDePoro-l));
X insertion du nouveau phenotype
Phen(i,NSTRAT+2)« Phen_possible(newPoroI,2);
X correspondance en binaire
newChrom_poro(1,:) = traduction(2,newPoroI,l:PRECI_poro);

ChromCi,((NSTRAT*PRECI_matID+ l*PRECI_poro+l):...
(NSTRAT*PRECI_matID+ 2*PRECI_poro)))...
=newChrom_poro(l,:); X porosité maximale
end X condition d'un indice trop eleve pour un Matériau
if (Phen(i,NSTRAT+2) == porosite_max )
paschanger=0;
else
paschanger=l;
end X if (Phen(i,NSTRAT+2) == porosite_max )
end X while pas changer
end X if NSTRA"=1
end X boucle en i sur les individus

X modifier les matériaux

for i ■ l:m X boucle sur les individus
for j=l:NSTRAT X boucle sur les phenotype
paschanger=0;
while paschanger==0
X indice matériau trop grand
if (Phen(i.j) >= MatDernier) k (Phen(i,NSTRAT+j)~=porosite_max)
X nouveau phenotype aléatoire
 newMat=ceil(rand(l)*MatDernier);
'/, insertion du nouveau phénotype
Phen(i,j)=newMat;
'/, position dans la liste de traduction
position = find(Phen_possible==newMat);
'/, correspondance en binaire
newChromd, : ) = traduction(l,posltion(l),1:PRECI_matID);
'/, porosité maximale
Chrom(i,( ( (j-l)*PRECI_matID+l ):j*PRECI_matID ))=ne«Chrom(1,:);
[im, in]=size(Chrom) ;
[mi,ni]= size(Chrom);
if ni~=nc
errorOcorrectegene: dimension de Chrom')
end
end 70 condition d'un indice trop élève pour un Matériau
if (Phen(i,j) >= MatDernier) & (Phen(i,NSTRAT+])"=porosite_max)
paschanger=0;
else
paschanger=l;
end
end '/. while
end "/, boucle sur les phénétiques
end '/. boucle en i sur les individus

'/.Fin

B.3.7 Tronquer la séquence

'/, troque.m - réduit la taille du stratifié dans certaines conditions

'/.
'/, Cette fonction vérifie si l'individu possède une porosité maximale
"/, parmi ses gènes, si tel est le cas, le stratifié est alors
'/, tronque a la première valeur de porosité maximale. Les porosités
% des gènes subséquents sont alors remplacés par une porosité maximale.
y.
% Syntaxe: [Chrom] = tronque(Chrom,PRECI_matID,PRECI_poro,NSTRAT)
y.
'/, Paramètre d'entrée :
y.
"/. Chrom - Matrice contenant la description des gènes
'/. sous la forme de chromosomes chaque ligne
*/, Correspondant a un individu. Le nombre de
■/, précision pour chaque gène est indique
'/. par PRECI_matID et PRECI_poro
'/.
7. Phen - table contenant les phénotypes des individus
'/. chaque ligne correspond a un individu
'/, chaque colonne correspond à une variable
'/. la première moitié des variables correspondantes aux
'/, matériaux qui composent les couches,
'/, la seconde moitié correspond a la porosité
'/. de chacune des couches
'/.
'/, Phen Matrice contenant la description des individus
7. Sous la forme de décimale.
7.
7. PRECIjnatlD Précision en nombre de bit pour le type de matéri
•/.
*/. PRECI_poro Précision en nombre de bit pour la porosité
'/.
'/. porosite_max porosité maximale admise après quoi
y. les porosités supérieures sont considéré
7. comme une couche vide
7,
'/. NSTRAT Nombre de couches dans le stratifié
7,
'/. MatDernier dernier matériau possible dans la liste avec le
•/. matériau vide
7,
'/, traduction Tableau contait la correspondance entre la
y. Représentation binaire et réelle.
7.
'/, Phen_possible tableau des différents phénotypes possible
7.
/.
'/, Paramètre de sortie:
y.
'/. Chrom Matrice contenant la description des gènes
y. sous la forme de chromosomes
7 chaque ligne correspondante a un individu,
y. Une fois, tronque.
y.
% Phen Matrice contenant la description des gènes
7, sous la forme de phénotypes
7 chaque ligne correspondant a
y. un individu, une fois tronque.
y.
'/. Auteur : Charles Villemure
'/, Historique: 31.08.05 : création de la fonction.

y,***** *************************************************************
function [Chrom,Phen] = tronque(Chrom,Phen,PRECI_matID,PRECI_poro,..
porosite.max,NSTRAT,MatDernier,...
traduction,Phen_possible)
y„******************************************************************.
[NIND,n]=size(Phen) ; */, n = nb de variables
7:1x1:1x1x1x1x1x1x1. b i n a i r e 7:1:1x1:1:1:1:1:1:1:1:1x1:1:1:1.

'/. P o r o s i t é maximale en format b i n a i r e

Poro_max=zeros(l,PRECI_poro);
Poro_max(l,1)=1;

'/, Porosité minimale en format binaire

Poro_min=zeros(l,PRECI_poro);

'/. Matériaux vide dans la liste ( dernier )

Mat_vide=zeros(l,PRECIjnatlD);
Mat_vide(l,l)=l;

vmxaitvavavmtravavavxanuivava'a
for i=l:NIND
for j=l:NSTRAT
% Si la porosité max est atteinte
if Phen(i,NSTRAT+j)== porosite_max;
X Boucle sur les couches subséquentes

'/, remplacer les couches subséquentes par des couches vides

for k=j:NSTRAT
'/, derniers matériaux dans la liste tableDmat — > couche vide
Phen(i,k)=2~PRECI_matID;
X imposition porosité maximale
Phen(i,NSTRAT+k)=porosite_max;
'/, couche vide — > dernier matériaux
Chrom(i,(PRECI_matID*(k-l)+l):(PRECI_matID*(k)))=Mat_vide;
7, porosité maximale
ChromCi,(NSTRAT*PRECI_matID+PRECI_poro*(k-l)+l):...
(NSTRAT*PRECI_matID+PRECI_poro*(k)))=Poro_max;
end X boucle sur les couches subséquentes
end
end '/, boucles sur les couches du stratifié NSTRAT
end 7, boucle sur individus : INDS
"/. Fin

B.3.8 Résolution thermofluide

y. naturelle.m - code CFD

'/. cette fonction calcule la distribution de température
'/. ainsi que la distribution de vitesse pour un stratifié
'/. donné en fonction des valeurs de porosités
'/, et le type de matériau de chacune des couches.
X
*/. Syntaxe: [ObjV,Dbase,T,v,Prop,strat,algo]
'/. = naturelle (strat, algoPhen, NSTRAT, Ra, Dbase, L, contraintes ,graph,
y.
% Paramètre d'entrée:
y. strat - propriétés du stratifié opt.
y. algo - propriétés de l'AG
y. Phen - table contenant les phénotypes des individus
X chaque ligne correspond a un individu
y. chaque colonne correspond a une variable
7, la première moitié des variables
y. correspond aux matériaux qui composent les couches
X la seconde moitié correspond
X a la porosité de chacune des couches
y. NSTRAT - Nombre de couches dans le stratifié
7. Ra - Nombre de Rayleigh : moteur de l'écoulement
7, Dbase - Base de donnée des calculs antérieurs
7. allure de la base de données
7 ( une ligne de la basse de données)
 129

"/, Dbase=[matériaux des couches, porosité des couches,...

'/, Ra, L, masseMax, coutMax, Ob j V, Tmax, ]
7. L - Longeur de l'échangeur L/H
7. contraintes - valeurs des contraintes
'/. contraintes de [masseMax, coutMax]
"/. graph - Paramètre optionnel :
5S indicateur pour tracer
'/, ou non les graphiques de température et vitesse
'/, si graph == 0 -> aucun graphique trace
'/. si graph == 1 -> graphiques traces
y.
'/, paramètre de sortie:
'/. strat - propriétés du stratifié opt.
*/. algo - propriétés de l'AG
'/. OptVal - valeur de référence sur laquelle repose l'optimisation du
'/. design. Celle-ci correspond a la valeur maximale de
% température atteinte dans la distribution de température à
°/. laquelle deux pénalités pénalités sont, ajoute, afin de
*/. répondre a des contraintes de : coût et de poids du
'/, Stratifié.
'/, T - Utilise pour un seul individu ( Phen == vecteur )
°/, Distribution de température de l'individu correspondant
'/. v - Utilise pour un seul individu ( Phen == vecteur )
'/. Distribution de vitesse de l'individu correspondant
'/. Prop - Propriété des individus:Prop = [masse_inds,cout_inds,Tmax]

y. Auteur : Charles Villemure

*/. 15.12.2006
y.
7,**************************************************************************
function [ObjV,Dbase,T,v,Prop,strat,algo] = naturelle(...
strat,algo,Phen,...
NSTRAT,...
Ra,Dbase,L,...
contraintes,graph,visu)
7a********************************************************* *****************
% Vérifier les paramètres d'entrer
if nargin == 8
graph=0;
visu=0;
elseif nargin »■ 9
visu=0;
elseif nargin < 8
errorCll manque des entrées a la fonction');
ond

ymxmxmi7xan%y:avhyxatvxauvxay:ayxavxayxi.
'/, Paramètres variable de la fonction et du maillage '/,
y^vxa^xayxayxayxixayxayxayxyxayxayxavx/x/xi.
'/. NY est fonction de Ra
°/e*******************************************
if Ra <= lel4
NY = 30; "/. Discrétisation en "y"
elseif Ra > lel4 & Ra <= lel5
NY = 60;
 130

elseif Ra > lel5 k Ra <= lel6

NY = 120;
elseif Ra > lel6 k Ra <= lel7
NY = 240;
elseif Ra > lel7 & Ra <= lel8
NY = 480;
end

'/, maillage et dimensions du domaine

H-l;
%L = 0.5; '/. Épaisseur de l'agencement (Épaisseur REELLE)
ep = ceil(NY*L/NSTRAT) ; '/, servant a multiplier le nombre
*/, de sous-couches par couche réelle
if ep < 5 '/, nombre de sous-couches minimales
ep = 5;
end
NLAYERS = ep*NSTRAT; '/, Discrétisation en "x"
'/. : nombre de couches (pour L'ALGORITHME)
DoverH = 0.001; '/, Diamètre des pores
'/. addimentionnel.
delta_y = 1/(NY-1); '/, Incrément en y
delta_y_2 = delta_y/2; '/. Incrément en y (moitié)
dx - L/(NLAYERS-1); % Incrémente en x

% Paramètre de convergence et fonction objectif

itermax=100; '/, nombre d'itération,

'/. maximal pour atteindre
w_T = 0.55; % facteur de relaxation pour,
% la température ( valeur optimsee )
w_v = 0.55; '/, facteur de relaxation pour la vitesse,
'/. ( valeur optimsee )
tolerance=10e-4; '/, tolérance sur la convergence,
'/, de la distribution,
Z de Température et de vitesse
porosité_max=l.0 ; '/, porosité maximale admise 100
Z correspond a une couche vide
massera = contraintes (1); '/. Masse que l'agencement
'/. ne doit pas dépasser
coutMax = contraintes (2); '/, Coût que l'agencement
'/, ne doit pas dépasser
TMaxlimite ■ contraintes(3); y, température limite
% ne doit pas dépasser
alpha=contraintes(4); y, paramètre de la fonction de pénalité
% applicable seulement sur TMaxlimite

coefl = 0.0; '/, Coefficient de pénalité de la mass

coef2 = 0.0; '/, Coefficient de pénalité du coût

'/. initialisation des variables

i=0;
j=0;
 131

n=0;

keq=zeros(l .NLAYERS) ; '/, conductivité équivalente au point P

keqE=zeros(l,NLAYERS) ; '/, conductivité équivalente au point Est
keqW=zeros(l,NLAYERS) ; '/. conductivité équivalente au point West

T = zéros(NY,NLAYERS) ; "/, distribution de Température dans le Composit

T_ancien=zeros(NY,NLAYERS) ; '/. distribution de Température
'/. dans le Composit itération précédente
v = zéros(1,NLAYERS); % vitesse a l'itération précédente
v_ancien=zeros (1, NLAYERS) ; '/. vitesse pour une couche donnée
int_T=zeros(l,NLAYERS) ; '/. intégrale de la vitesse

aP= zéros (1, NLAYERS) ; '/. coefficient de TP en Volume de contrôle.

aN= zéros ( 1, NL AYERS ) ; '/. coefficient de TN en Volume de contrôle.
aS= zerosO .NLAYERS) : "/, coefficient de TS en Volume de contrôle.
aW= zéros (1, NLAYERS) ; '/. coefficient de TW en Volume de contrôle.
aE= zéros(1,NLAYERS) ; '/. coefficient de TE en Volume de contrôle.
r- zéros(1 .NLAYERS) ; '/. coefficient de TE en Volume de contrôle.
A= zéros (NLAYERS, NLAYERS); '/. A*T = r
aire_Ty = zéros (NY,NLAYERS) ; '/. aire sous la courbe pour un delta y
int_Ty = zéros(NLAYERS); % intégrale de l'aire sous la courbe

%*************************************************************************
iterHist=l;

°/ 0 * ************************************************************************
yx/x/x/x/x/x/x/x/x/x/x/x/x/x/x/xa
'/, Propriété du stratifié poreux "/,
yxixixixixmxixixixixixixixixixix
y,*************************************************************************
'/, Déterminer les dimensions de la population
[NIND,NVAR]=size(Phen) ;

'/. Propriétés du Stratifié: Porosité, Densité, Conductivité, Coût

[densité, k,porosité, cout_inds, masse_inds] ...
= propriété(Phen,NSTRAT,ep,dx,porosite_max);

'/. Sortie de la fonction Prop » [masse,coût] pour chaque individu

Prop(:,1)=masse_inds;
Prop(:,2)=cout_inds;

coutMax_inds=ones(NIND,1).*coutMax;
masseMax_inds=ones(NIND,l).*masseMax;

vx/x/x/x/x/x/x/x/x/x/x/x/x/x/x/xa
"/, Propriété du Stratifié poreux %
yx/x/x/x/x/x/x/x/x/x/x/x/x/XÂt'/xa
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
[ObjVhist,Tmaxhist]"historique(0,Phen,Dbase,NSTRAT,Ra,L.masseMax,coutMax);

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y. Calculs des conductivités intermédiaires '/,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y**********************************************************************

'/, Déterminer les conductivités intermédiaires (sous-couches)

[keq.keqE.keqW] = conductivité (k,NIND,NLAYERS,porosité);
y.[keq,keqE,keqW] = conductivité_bejan(k,NIND,NLAYERS,porosité) ;

% Déterminer les perméabilités (sous-couches)

[permea] = perméabilité (porosité,DoverH,NIND,NLAYERS);

avant=0;
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
X Resoulution de la Distribution de Température '/,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%**********************************************************************
X identifiant d'une porosité nulle
poro_nulle=0;

'/, Déterminer si une couche possède une porosité égale a 1

for inds=l:NIND X boucle sur les individus

X
'/. Verfier si l'individu a déjà été calculé auparavant
X

'/, condition : pas la première ligne de Dbase (-1)

'/, condition : pas de graphique demandé à naturelle

if (ObjVhist(inds) "- -1 & graph~=l) && isfinite(ObjVhist(inds))

'/. On récupère la valeur ObjV calcule app.

ObjVUnds.l) = ObjVhist(inds.l);

'/. ajouter aux propriétés du stratifié la température maximale

Prop(inds,3) = Tmaxhistdnds, 1) ; '/, récupère la valeur de Tmax

% mettre à jour le nb récupération dans la basse de donne

algo.nbdbase=algo.nbdbase+1 ;

else '/, on effecture le calcul itératif

yx/x/x/x/x/x/x/x/x/x/x/x/.yx/x/xm
XXXXX début de la résolution du profil de température XXXXXXXX
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

y. Vérifier si la porosité est totale

for j-i:NLAYERS
if ( porosité(inds.j) == porosite_max)
Nadia = j ;
break
else
Nadia = NLAYERS;
end
end
 X vecteur des dy utilise pour calculer l'air sur la courbe de Ty a x=0
for i= 2:NY-1
vect_dy(i,l)=delta_y;
end
vect_dy(1,1)=delta_y_2;
vect_dy(NY,1)=delta_y_2;

X Resolution de la température pour un individu X

X
'/, initialisation
iter=l;
conv=0;
T(l,:)=0.5; X afin d'évité une vitesse nulle a la première itération

while (iter < itermax)

X Critère qui provoque l'arrêt lorsque la convergence est atteinte

% avant que le nombre maximum d'itération soit effectue

if ( conv == 1 ) 7,si convergence atteinte conv = 0, sinon conv = 1 ;

if iter==itermax
warningCil manque de latitude pour la convergence en T ou v ' ) ;
break;
end
break;
end
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
'/. intégration du profil de température tVXIXIXIXIXIXIXIXIXI'l,
xxxxxxxxxxxxxxxxxmxmxxxxxmxxxxxxxxxxxxxxxxxxxxxxxxmxxxxmmxxx
X for ( j=l:NLAYERS )
7. aire_Ty(: ,j)=T(: ,j) .*vect_dy;
X int_Ty(j)=sum(aire_Ty(:,j));
% end
X
X X vitesse en "x"
X for ( j=l:NLAYERS )
X Xv(j)=permea(inds,j)*Ra*T(l,j); 7. modèle simple
X v(j)=permea(inds, j)*Ra*int_Ty(j) ; 7. modèle plus complexe
X Xv(j)=1000; 7. convection forcée
7. end
X
y.
int_Ty2(j)=int_Ty(j);
for j=l:NLAYERS
7int_Ty(j) = trapèze(transpose(T(:,j)),transpose(vect_dy));
int_Ty(j) = trapezeC(transpose(T(;,j)),transpose(vect_dy));

v(j)=permea(inds,j)*Ra*int_Ty(j);
end X boucle sur les couche NLAYERS
X
X résolution du profil de température
X

X Boucle sur nombre d'intervalles en Y en partant du haut

for(i=NY:-l: 1) X Boucle sur nombre d'intervalles en Y
 /, NY-1 car la première ligne est fixe à T=0 — > CL.

if i==l '/, demi-volume de contrSle a l'extrémité

dy=delta_y_2;
else
dy=delta_y;
end

'/, intérieur du domaine

for(j=2:(NLAYERS-l))

% cas au centre
aP(j) = -dy/dx*keqE(inds,j) -dy/dx*keqW(inds,j) -dx*v(j);
aN(j) = 0;
aS(j) = dx*v(j);
aW(j) = dy/dx*keqW(inds,j);
aE(j) = dy/dx*keqE(inds,j);

if i==(NY)
r(j) = 0 ; % condition limite entre
else
r(j) = -aS(j)*T(i+l, j) ; '/. partout ailleurs
end
end */, % boucle en j — > largeur du stratifié

y. condition limite x = 0; flux constant : dT/dx=cste

J"li
aP(j) = -dy/dx*keqE(inds,j) -dx/2*v(j);
aN(j) = 0;
aS(j) - dx/2*v(j);
aVKj) = 0;
aE(j) = dy/dx*keqE(inds,j);

if i==(NY)
r(j) = -dy ;
else
r(j) = -dy -aS(j)*T(i+i,j)j
end

7. condition limite x = L/l; flux constant : dT/dx=null

j=NLAYERS;
aP(j) = -dy/dx*keqW(inds,j) - dx/2*v(j);
aN(j) = 0;
aS(j) - dx/2*v(j);
aW(j) = dy/dx*keqW(inds,j);
aE(j) = 0;

if i==(NY)
r(j) = 0;
else
r(j) = -aS(j)*T(i+l,j);
end

"/, condition pour une porosité totale

*/, paroi adiabatique sur la couche Nadia
 if Nadia"=NLAYERS
j=(Nadia-l);
aP(j) = -dy/dx*keqW(inds,j) -dx/2*v(j);
aN(j) = 0;
aS(j) = dx/2*v(j);
aW(j) = dy/dx*keqW(inds,j);
aE(j) = 0;

if i==(NY)
r(j) = 0;
else
r(j) = -aS (j)*T(i+1,j)
end

f or(j=Nadia:NLAYERS)
aP(j)=i;
aE(j)=0;
aW(j)=0;

r(j)=0:
(MU]

end % Nadia "= NLAYERS

% résolution de la matrice tridiagonale :

T_ligne=tridiagC(aW(l:NLAYERS),aP,aE(l:(NLAYERS-i)),r);
V.T_ligne=tridag(aW( 1 : NLAYERS), aP, aE( 1 : (NLAYERS-1) ) , r) ;

'/, position le vecteur T dans sa direction physique

T(i,:)=transpose(T_ligne);

end "/, boucle en i — > hauteur du composit

'/, Test de la convergence : Résidus %

'/, sommation du carre des résidus

il iter <=5
conv_residus(iter)=sum(sum(abs(résidus)))/sum(sum(abs(Rap_normal)));
max_conv_residus=max(conv_residus);
else
conv_residus(iter)=sum(sum(abs(résidus)))...
/sum(sum(abs(Rap_normal)))/max_conv_residus;
end
il iter >6
îor j i=l : 5
conv_residus(ji)=conv_residus(6);
end
end

'/, convergence de la vitesse

 R_vitesse_brut(iter)=sqrt(sum((v-v_ancien)."2));

if iter < 5
conv_v=0;
R_normal_vitesse=max(R_vitesse_brut);
R_vitesse(iter)=R_vitesse_brut(iter)/R_normal_vitesse;
else
R_vitesse(iter)=R_vitesse_brut(iter)/R_normal_vitesse;
if R_vitesse(iter) < tolérance
conv_v=l;
else
conv_v=0;
end
end

'/.
'/, convergence de la Température

R_tempera_brut(iter)=sqrt(sum(sum((T-T_ancien)."2)));

if iter < 5
conv_T=0;
R_normal_t empera=max(R_t empera_brut);
R_tempera(iter)=R_tempera_brut(iter)/R_normal_tempera;
else
R_tempera(iter)=R_tempera_brut(iter)/R_normal_tempera;
if R_tempera(iter) < tolérance
conv_T=l;
else
conv_T=0;
end
end

if conv_T ==1 kk conv_v==l

conv=l;
else
y.
'/, relaxation
T=T_ancien + w_T.*(T-T_ancien);
v=v_ancien + ¥_v.*(v-v_ancien);
end

'/, mise a jour des distributions

v_ancien=v;
T_ancien=T;

% température maximale dans le stratifié

T_max(iter)=max(max(T));

%
'/, Visualisation a l'écran :
'/, graph de diff_T, diff_v et T_max en fonction des itérations.
% %
if visu==l
figure(inds)
 subplot(2,2,1);plot(loglO(T_max),'ro'); xlabeK'température maximale');
ylabel('loglO');
text(0.5,0.90,['T_max = ', num2str(T_max(iter))],'Units','normalized');
drawnow;

subplot(2,2,2);plot(loglO(R_tempera),'ro');...
xlabeK'résidus température conserved');
ylabel('loglO');
drawnow;

subplot(2,2,3);plot(loglO(conv_residus),'ro');
xlabeK'résidus température');
ylabel('loglO');
drawnow;

subplot(2,2,4);plot(loglO(R_vitesse),'ro');
xlabeK'différence de vitesse');
ylabel('loglO');
drawnow;
end
7. 7.
t

iter=iter+l ; 7. intégrateur
end 7, While pas convergence
7. fin de la Resolution de distribution de Temperaure

utixavata'/xa
VI. Optimisation 7.

ininnvxauu
7.display(iter) ;

7. Valeurs a minimiser
Tmax(inds,l) = max(max(T));
ObjVavantCinds, 1) = Tmax(inds,l) ;

7. ajouter aux propriétés du stratifié la température maximale

Prop(inds,3) = Tmax(inds,l);
Temp=Tmax(inds,1);

7. Propriété du stratifié pour l'individu "inds"

UNmasse = masse_inds(inds,l);
UNcout = cout_inds(inds,l);

7. Contraintes imposées a l'ensemble des individus

7. Pénalise l'individu si les caractéristiques excèdent les contraintes
if TMaxlimite==0.0
[ObjVdnds, 1) ,strat,algo]=penalite(strat,algo,ObjVavant(inds, 1), .. .
UNmasse,UNcout,...
masseMax.coutMax);
else
[ObjVdnds, 1) > strat,algo]=penalitealpha(strat,algo,ObjVavant(inds,l),. . .
UNmasse,UNcout,Temp,...
masseMax,coutMax,TMaxlimite,alpha);
 end '/. if TMaxlimite==0.0

V, Déterminer les dimension de la base de donne

[tous,Ncaract]=size(Dbase);

'/, Ajout a la base de donne

Dbase(tous+l,:)=[Phen(inds,:),Ra,L,masseMax,...
coutMax,ObjV(inds,l),Tmax(inds,1)];
'/, mettre à jour le nb d'évaluation
algo.eval=algo.eval+1;
end '/, if de histoire ObjVhist
°/B ************************************************************************
iv/:a°avava
% Graphique '/,
mv/x/xava
%* ************************************************************************
if (graph==l)

for i=l:NY
hauteur(i)=(i-l)*(l/(NY-l));
end

for i=l:NLAYERS
largeur(i)=(i-l)*(l/(NLAYERS-l))*L;
pos_x(i)=(i-l)*L/(NLAYERS-l);
end

for i=l:(itermax-1)
nb_iter(i)=i;
end

figure (20)
meshdargeur,hauteur,T);

figure (21)
plot(pos_x,v);

figure (22)
contourdargeur,hauteur,T) ;

end '/, graph ==1

%************************************************************************
end '/, résolution pour les individus
% Fin

B.3.9 Conductivité équivalente

'/, conductivité.m - calcule la conductivité équivalente du milieu

'/, Cette fonction calcule la conductivité équivalente du milieu.
'/, Cette conductivité est déterminée selon une répartition harmonique.
y.
'/. Syntaxe:
'/. [tableDdensite,tableDk,tableDmat.tableDcout,nbDmat]= matériaux (CHX)
'/.
'/, Paramètre d'entrée:
'/.
'/, k - Conductivité thermique pour chacune des sous-couches
'/. composant le stratifiée
'/. NIND - nombre d'individus qui compose une population
'/, NLAYERS - nombre de sous-couches qui composent le stratifié
'/, porosité - vectoeur comprenant la valeur de la porosité pour
% chacune des sous-couches
% Paramètre de sortie:
y.
% keq - Conductivité thermique équivalente au point central
'/ keqE - Conductivité thermique équivalente au point central
'/. keqW - Conductivité thermique équivalente au point central
y.
y.
'/o Auteur : Charles Villemure
•/. 20.12.2006

function [keq,keqE,keqW] = conductivite(k,NIND,NLAYERS,porosité)

y^********* ************************************* ****************************
keq=zeros(NIND,NLAYERS);
keqE=zeros(NIND,NLAYERS);
keqW=zeros(NIND,NLAYERS) ;

'/. Point P : centre : calcul de la conductivité équivalente

for i=l:NIND % boucle sur les individus
for j=l:NLAYERS '/, boucle sur les couches
keq(i,j)=k(i,j)*((l-porosite(i,j))/(l+porosite(i,j)));
end X boucle sur les couches
end '/, boucle sur les individus

% Point : Est
for i=l:NIND "/, boucle sur les individus
for j=l:NLAYERS-l "/, boucle sur les couches

°/, cas ou les deux conductivités équivalentes sont nulles

if (keq(l.j) == 0 kk keq(i,j+l) == 0)
keqE(i,j)=0;
else
'/, moyenne harmonique
keqE(i,j)=(2*keq(i,j)*keq(i,j+l))/(keq(i,j)+keq(i,j+l)) ;
end
end V, boucle sur les couches
keqE(i,NLAYERS)=keq(i,NLAYERS);
end y. boucle sur les individus

y. Point : West
for i=l:NIND "/. boucle sur les individus
for j=2:NLAYERS % boucle sur les couches
 '/, cas ou les deux conductivités équivalentes sont nulles
if (keq(i.j) == 0 kk keq(i,j-l) == 0)
keqW(i,j)=0;
else
keqW(i,j)-(2*keq(i,j)*keq(i,j-l))/(keq(i,j)+keq(i,j-l));
end
end '/, boucle sur les couches
keqW(i,l)=keq(i,l);
end '/. boucle sur les individu
t Fin

B.3.10 Perméabilité

'/, perméabilité.m - Calcul de perméabilité du milieu poreux

7.
"/, Cette fonction calcule des conductivités intermédiaires
'/, pour les points composant les sous-couches.
'/, ces conductivités sont déterminer selon % une répartition harmonique.
y.
'/, Syntaxe: [permea] = permeabilite(porosite,DoverH,NIND,NLAYERS)
y.
'/. Paramètre d'entrée:
y.
'/, porosité - vecteur comprenant la valeur de la porosité pour
'/, chacune des sous-couches
y.
y, DOoverH - diamètre des pores (addimentionel) D/H
y.
'/, NIND - nombre d'individus qui composent une population
y.
V, NLAYERS - nombre de sous-couches qui composent le stratifié
y.
'/, Paramètre de sortie:
y.
'/, permea - perméabilité pour chacune des sous-couches
y.
'/, Auteur : Charles Villemure
Z 20.12.2006
x

function [permea] - permeabilite(porosite,DoverH,NIND,NLAYERS)

y*********************************************************************

'/. calcul des perméabilités à partir des porosités (forme adimensionnelle)

for i=l:NIND X boucle sur les individus
for(j=l:NLAYERS)
permea(i,j) = porosited, j)*DoverH*DoverH*l/32;
end y,f or (j-1: NLAYERS)
end '/. boucle sur les individus
X Fin

B.3.11 Résolution d'un système triagonal

X tridiag.m - cette effectue la résolution d'un système matriciel

X tridiagonal.
X
'/, Syntaxe: u=tridag(a,b,c,r)
X
% Paramètre d'entrée:
X a - vecteur contenant les valeurs de diagonale inférieure
'/.
X b - vecteur contenant les valeurs de diagonal central '/,
X c - vecteur contenant les valeurs de diagonal supérieur X
X r - vecteur contenant les valeurs des termes source
X
X paramètre de sortie:
X
X u - vecteur contenant la solution du système A*U=r;
X
X Système :
X
X bl cl 0 ... I
X a2 b2 c2 0 I
X A- 0 a3 b3 c3 0 I A*U=r;
X I
X 0 0 aN bN I
X
X Auteur : Charles Villemure
X Historique : 27.01.06 : création de la fonction.
X Référence : http://wuw.princeton.edu/~lam/course/tridag.rn
X
function u=tridag(a,b,c,r);
X

if nargin<4
disp('nécessite a, b, c et r en entrée.')
end
n=length(b);
if length(a)==n-l
a» [0 a] ;
elseif length(a)==n
a(l)«0;
end
if length(b)~=length(r)
dispCles vecteurs doivent avoir la même dimension ')
end
if length(c)<n-l
dispCLe vecteur c est trop court')
end
'/. Numerical Recipes, Press, Flannery, Teukolsky and Vetterling,
'/, Cambridge Univerisity Press (1986), section 2.6, page 40.
'/. Simple LU décomposition algorithm.
bet=b(l);
if bet==0
dispCbl ne peut être égale à 0 . . . ' ) ;
bet=eps;
end
u(l)=r(l)/bet;
for j=2:n
gam(j)=c(j-l)/bet;
bet=b(j)-a(j)*gam(j);
if bet==0
'/.dispC... probablement, baisions de pivoter ' ) ;
bet=eps;
end
u(j)=(r(j)-a(j)*u(j-l))/bet;
end
for j=n-l:-l:l
u(j)=u(j)-gam(j+l)*u(j + l) ;
end

B.3.12 Intégration d'une fonction

/**"/. trapèze.m - Cette fonction calcule la l'intégral

1, d'une fonction numériquement par la méthode des trapèzes.
%
% Syntaxe: [intT] = t r a p è z e ( T , y )
Z
'/. Paramètre d ' e n t r é e :
'/, T - vecteur contenant les valeurs numériques a intégré
y.
"/, y - axe par rapport a lequel la T sera intègre
y.
'/. paramètre de sortie:
y.
'/, intT - valeur de l'intergrale
y.

% Auteur : Charles Villemure

X Historique: 27.01.06 : création de la fonction.**/

#include "mex.h"
void mexFunctiondnt nbArgsOut, mxArray *sortie[],
int nbArgsIn, const mxArray «entrée [])
{ /** début fx **/

int longT;
int un;
int i;
double *inT;
double *T;
double *yt;

un=l;

//inT = mxGetPr(raxCreateDoubleArray(l,un,mxREAL));

/ * * Longeur de l ' e n t r é e 1, T * * /
longT = mxGetN(entree[0]);

/** Créer la matrice de sortie **/

s o r t i e [0] = mxCreateDoubleMatrix(l,un,mxREAL);
inT = m x G e t P r ( s o r t i e [ 0 ] ) ;

/ * * a l l e r chercher l ' a d r e s s e des e n t r é e s * * /

T = mxGetPr(entrée[0]);
yt » m x G e t P r ( e n t r e e [ l ] ) ;

*inT=0;

for(i=l;i<longT;i++ )
{
*inT=*inT+(T[i-l]+T[i])/2*yt[i-l] ;
}
}/**fin fx**/

B.3.13 Pénalité de la fonction objectif

"L pénalité.m - pénalise la fonction objectif

X
V, Cette fonction modifie la valeur de l'objectif
'/. en lui attribuant une pénalité selon le critère choisi
% et les valeurs de contrainte de masse et de coût.
y.
% Syntaxe: [ObjVmod,strat,algo]
Z = penalite(strat,algo,ObjV,masse,coût,masseMax,coutMax)
y.
'/, Paramètre d'entrée:
'/. strat - propriétés du stratifié opt.
'/o algo - propriétés de l'AG
'/, ObjV - Valeur a optimisé : Tmax
'/, caractéristique - Caractéristique du stratifié : [masse,coût]
y. contraintes - Contraintes de l'optimisation : [masseMax,coutMax]
y, si masseMax == 0
'/. -> aucune optimisation sur la masse
'/, si coutMax == 0
"/, -> aucune optimisation sur le coût
"/, Paramètre de sortie :
7, ObjVmod - valeur a optimiser modifie
'/.
'/, Auteur : Charles Villemure
•/. 20.12.2006

y,* ************************************************************** ***********

function [ObjVmod,strat,algo]=penalite(strat,algo,...
ObjV,masse,coût,...
masseMax,coutMax)
7,* *************************************************************************
7. Coefficient devant la parabole pour la pénalité
coeff=20;
% ratio des carac / cont
if masseMax ~=0
ratio_masse = masse / masseMax;
else
ratio_masse =0;
end

if coutMax ~=0
ratio_cout = coût / coutMax;
else
ratio_cout=0;
end

if r a t i o _ c o u t < 0 I r a t i o n n a s s e < 0
error('erreur dans la pénalité au niveau des ratios');
end

'/, si la masse excède la masse maximale

if rationnasse > 1 k masseMax "■ 0
contr_masse = coeff*(ratio_masse-l)"2 + 1;
7, valeur a multiplier
elseif ratio_masse < 1 | masseMax == 0
contr_masse = 1;
end '/, if r a t i o n n a s s e

'/, si le coût excède le coût maximal

if ratio_cout > 1 k coutMax "= 0
contr_cout = coeff*(ratio_cout-l)~2 + 1;
elseif ratio_cout < 1 | coutMax == 0
contr_cout = 1 ;
end °/, if ratio_masse

'/. Pénalisation de ObjV

ObjVmod = ObjV * contr_masse * contr_cout;

'/. fin
B.3.14 Classement

'/, ranking.m - trie les individus selon leur performance

•/.
7. Cette fonction trie les individus selon leur degré de performance
7. relative au sein de la population. ( minimisation de la fonction
*/. objectif)
•/.
'/. Syntaxe: FitnV = rankingCObjV, RFun, SUBPDP)
7.
y,
7. Paramètres d'entrée :
7. ObjV - Vecteur colonne contenant la valeur de la fonction objectif
7, de chaque individu de la population.
'/. RFun - Paramètre optionel. Si RFun est un scalaire, un trie
'/, linéaire est réalisé et le scalaire indique la pression de
7. sélection.
7. SUBPOP - Paramètre optionnel. Nombre de sous-populations
•/.
7. Paramètres de sortie:
'/, FitnV - Vecteur colonne contenant la valeur associocié à la
7. Performance relative d'un individu.
7.
7, Auteur: Hartmut Pohlheim,Carlos Fonseca
7. Adaptation: Charles Villemure
7. 15.20.2006
7.
7**************************************************************************
function FitnV = rankingCObjV, RFun, SUBPOP);
y 0 * * * * * * * * * * * * * *************************************************************
7. Déterminer l a t a i l l e de l a population
[Nind,ans] = size(ObjV);

7. Vérifier les paramètres d'entrée

if nargin < 2, RFun = [] ; end
if nargin > 1, if isnan(RFun), RFun = []; end, end
if prod(size(RFun)) == 2,
if RFun(2) == 1, NonLin = 1;
elseif RFun(2) == 0, NonLin = 0;
else errorCLe choix du type de méthode indiqué par 0 ou 1'); end
RFun = RFun(l);
if isnan(RFun), RFun = 2; end
elseif prod(size(RFun)) > 2,
if prod(size(RFun)) ~= Nind
errorCLes longueurs de ObjV et RFun ne concordent pas');
end
end

if nargin < 3, SUBPOP = 1; end

if nargin > 2,
if isempty(SUBPOP), SUBPOP = 1;
elseif isnan(SUBPOP), SUBPOP - 1;
elseif length(SUBPOP) '- 1
errorCSUBPOP doit être un scalaire');
 ont!
end

if (Nind/SUBPOP) "- fix(Nind/SUBPOP)

e r r o r C O b j V et SUBPDP ne concordent p a s ' ) ;
end
% détermine le nombre d'individus par sous-population
Nind = Nind/SUBPOP;

'/. Identifie la fonction de trie

if isempty(RFun),
"/, trie linéaire avec une pression de sélection de 2 ( défaut)
RFun = 2*[0:Nind-l]'/(Nind-l);
elseif prod(size(RFun)) == 1
if NonLin == 1,
"/, trie non-linéaire
if RFun(l) < 1
e r r o r C L a pression de sélection doit être sup. à 1 ' ) ;
elseif RFun(l) > Nind-2
errorC'trop grande valeur de la SP ' ) ;
end
Rootl = roots([RFun(l)-Nind [RFun(l)*ones(l,Nind-l)]]);
RFun = (abs(RootKD) * ones(Nind.D) .* [(0:Nind-l)'] ;
RFun = RFun / sum(RFun) * Nind;
else
'/. linear ranking with SP between 1 and 2
if (RFun(l) < 1 I RFun(l) > 2 ) ,
errorCEn trie linéaire, la SP doit être entre comprise entre[0,1] ' ) ;
end
RFun = 2-RFun + 2*(RFun-l)* [0:Nind-l] V(Nind-l) ;
end
end;

FitnV = [] ;

'/, Boucle sur les sous-pop

for irun = 1:SUBP0P,
V. copier les valeurs objectives de la population initiale
ObjVSub = ObjV((irun-l)*Nind+l:irun*Nind);.
'/. Déterminer s'il y a des NaN dans les valeurs de Fxobjectif
NaNix = isnan(ObjVSub);
Validix = find("NaNix);
% Trier seulement les valeurs numériques ( on minimise ) .
[ans.ix] = sort(-DbjVSub(Validix));
'/, Construire un vecteur des Fx Objectif en supposant que les NaN sont
% pire que les valeurs réelles. ( Incluant les infinity )
ix = [find(NaNix) ; Validix(ix)];
°/, La version finale
Sorted = ObjVSub(ix);
'/, Assigner la valeur de trie selon RFun
i = 1;
FitnVSub = zeros(Nind,1);
for j - [find(Sorted(l:Nind-l) "= Sorted(2:Nind)); Nind]',
FitnVSub(i:j) = sum(RFun(i:j)) * ones(j-i+l,1) / (j-i+1);
i -j+i;
 end
7. Un vecteur de trie,
'/, associé a la position dans la matrice population
[ans.uix] = sort(ix);
FitnVSub = FitnVSub(uix);
% juxtaposer les vecteurs de tous les sous-pop.
FitnV = [FitnV; FitnVSub];
end
'/. fin

B.3.15 Recherche locale

locale.m - cette fonction effectue une recherche locale

Syntaxe: [PhenMOD] = correctgenelocale(Phenmeil,Chrommeil,FieldD,...

PRECI_matID,PRECI_poro,...
NSTRAT.MatDernier,porosite_max,...
Ra,Dbase,L,contraintes)
Paramètre d'entrée:3

Phennmeil Table contenant les meilleurs individus d'une

Population. Chaque ligne correspond à un
individu chaque colonne correspond à une
variable. La première moitié des variables
correspondent aux matériaux qui composent
les couches. La seconde moitié correspond
à la porosité de chacune des couches

Chrommeil Matrice contenant la description des gènes

des meilleurs individus d'une population
(associé à la matrice Phenmeil)
sous la forme de chromosomes, chaque ligne
Correspondant a un individu.

FieldD Caracatéristique de la représentation

binaire d'un individu
( une colonne pour chaque individu )

PRECI_matID Précision en nombre de bit

pour le type de matériel

PRECI_poro Précision en nombre de bits

pour la porosité

NSTRAT Nombre de couches dans le stratifié

MatDernier Indice du dernier matériau possible

Paramètre de sortie:

PhenMOD Phénotypes modifiés

% Auteur : Charles Villemure
'/. 20.12.2006
y0**************************************************************************
function [Phenmeil,Chrommeil,Dbase,stumieux,strat,algo] = ...
locale(Phenmeil,Chrommeil,FieldD,...
PRECI_matID,PRECI_poro,...
NSTRAT,MatDernier,porosite_max,...
Ra.Dbase.L,contraintes,stumieux,...
strat,algo)
y o * *************************************************************************

X construire la table de traduction

[traduction,Phen_possible] = enigma(FieldD);

"/, nombre de porosités possibles

NbDePoro=length(Phen_possible(:,2));

'/, dimension du tableau des phenotypes

[m,n]=size(Phenmeil);

'/, dimension du tableau des phenotypes (binaire)

[me,ne]=size(Chrommeil);

'/, copie des individus originaux

PhenMOD=Phenmeil;
ChromMOD=Chrommeil ;

'/. variables de porosité

porosite=zeros(m,NSTRAT);
porosite=Phenmeil(:.NSTRAT+1:2*NSTRAT);

*/, une recherche dans le voisinage au niveau de la porosité

'/, est réalisé pour chacun des inds
'/. si la recherche et fructueuse,
"/. l'individus des remplacés par le résultat.

for i = l:m '/, boucle sur les individus

'A nombre de gènes qui subissent un amélioration local
nbgenes=ceil(rand(l)*NSTRAT);

'/, vecteur contenant les positions aléatoires des gènes

ordregenes=randperm(NSTRAT);

"/, vecteur contenant les gènes qui subissent un amélioration locale

indicegenes=ordregenes(l:nbgenes);

'/, matrice contenant les améliorations locales ( phenotypes locaux )

Phenlocal=zeros(nbgenes*2,2*NSTRAT);
Chromlocal=zeros(nbgenes*2,nc);

'/, la même distribution de matériaux est conservée (2 copies )

for igene=l :nbgenes
PhenlocaKigene, : )=Phenmeil(i , : ) ;
PhenlocaKigene+nbgenes, : )=Phenmeil(i, : ) ;
 149

ChromlocaKigene, : )=Chrommeil(i , :) :
ChroralocaKigene+nbgenes, : )=Chrommeil(i , : ) ;
end

for j = lmbgenes '/. boucle sur les gènes porosités

lequelgene=indlcegenes(j) ; '/, indice du gène a modifié
saporosite=zeros(l .NSTRAT) ; '/, porosité pour l'inds.
saporosite=porosite(i,:);
7. trouver la position du gène dans la liste enigma
if lequelgene==l
7. s'il agit du premier gène de porosité — > collone 3
if NSTRAT~=1
positionliste=find(Phen_possible(:,3)==saporosite(lequelgene));
else
positionliste=find(Phen_possible(:,2)==saporosite(lequelgene));
end 7. if NSTRAT~=1
else
7. s'il agit du premier gène de porosité — > collone 2
positionliste=find(Phen_possible(:,2)==saporosite(lequelgene));
end 7. lequelgene==l

if isempty(positionliste)
error('positionliste=0, recherche dans enigma infructueuse');
end

7. modification dans le voisinage

if positionliste «■ 1
7. position du gène supérieur
positionsup=positionliste+l;
7. position du gène supérieur ( 2e ordre )
positioninf=positionliste+2;
elseif positionliste == 2"(PRECI_poro)
7. position du gène supérieur
positionsup=positionliste-l;
7. position du gène supérieur ( 2e ordre )
positioninf=positionliste-2;
else
7. position du gène supérieur
positionsup=positionliste+l;
7. position du gène inférieur
positioninf=positionliste-l;
end 7opositionliste ~= 1

7. trouver la position du gène dans la liste enigma

if lequelgene==l
7. s'il agit du premier gène de porosité — > collone 3
if NSTRAT~=1
genesup=Phen_possible(positionsup,3);
geneinf=Phen_possible(positioninf,3);
binairesup=traduction(3,positionsup,l:PRECI_poro);
binaireinf=traduction(3,positioninf,1 :PRECI_poro);
else
genesup=Phen_possible(positionsup,2);
 geneinf=Phen_possible(positioninf,2);
binairesup=traduction(2,positionsup,l:PRECI_poro);
binaireinf=traduction(2,positioninf,l:PRECI_poro);
end %if NSTRAT"=1
else
"/, s'il agit du premier gène de porosité — > collone 2
genesup=Phen_possible(positionsup,2);
geneinf=Phen_possible(positioninf,2);
binairesup=traduction(2,positionsup,l:PRECI_poro);
binaireinf=traduction(2,positioninf,1:PRECI_poro);
end '/, lequelgene==l

"/, copie des vecteurs porosités

porositesup=saporosite;
porositeinf=saporosite;

% remplacer le gènes par un du voisinage

porositesup(lequelgene)=genesup;
porositéinf(lequelgene)=geneinf;

'/. Mettre a jour la matrice de Phénotype locaux

Phenlocal (j , NSTRAT+1:2*NSTRAT )=porositesup;
Phenlocal (j+nbgenes , NSTRAT+1:2*NSTRAT )=porositeinf;
sMat=NSTRAT*PRECI_matID;
sPoro=(lequelgene-1)*PRECI_poro;
somme=sMat+sPoro;
ChromlocaKj , (somme+1):(somme+PRECI_poro) )=binairesup;
ChromlocaKj+nbgenes , (somme+1):(somme+PRECI_poro) )=binaireinf;

•/, tester la concordance

test=bs2rv(Chromlocal,FieldD);
marcheti=sum(sum(test-Phenlocal));
if marcheti"=0.0
errorCerreur de concordance entre Chromlocal et Phenlocal');

end
end Y, boucle en j sur les gènes porosités

[checkNind,checkNvar]=size(Phenlocal);
if checkNvar~=NSTRAT*2
save('ErreurCheckNind')
error('ErreurCheckNind')
ond

y. Correction des gènes

[Phenlocal.Chromlocal] « correctgene(Phenlocal,Chromlocal,FieldD,...
PRECI_matID,PRECI_poro,NSTRAT,...
MatDernier,porosite_max);

y. Tronquer le stratifié si l'on rencontre une couche a porosité maximale

[Chromlocal,Phenlocal] = tronque(Chromlocal,Phenlocal,...
PRECI_matID,PRECI_poro,...
porosite_max,NSTRAT,MatDernier,...
traduction,Phen_possible);
 7. Évaluation de la population issus de la recherche locale
[ObjVlocal,Dbase,T,v,Prop,strat.algo] = n a t u r e l l e ( s t r a t . a l g o , . . .
Phenlocal.NSTRAT,...
Ra,Dbase,L,contraintes);

'/, Évaluation de la population initiale

[ObjV,Dbase,T,v,Prop,strat.algo] = naturelle(strat,algo,...
PhenmeiKi, :) ,NSTRAT,Ra,Dbase,L,contraintes) ;

'/, comparer la population issue de la recherche locale avec l'individu

[Elbest,lequelbest] = min(ObjVlocal);

% Vérifier si l'amélioration locale est meilleure if Elbest < ObjV

7. l'individu est remplacé
PhenmeiKi, : )=Phenlocal(lequelbest, : ) ;
ChrommeiKi, : )=Chromlocal(lequelbest, : ) ;
ObjV(i)=ObjVlocal(lequelbest);
stumieux=stumieux+l;
7.display(stumieux)
end '/. condition Elbest
end '/, boucle en i sur les individus

7. une recherche dans le voisinage au niveau des matériaux

1, est réalisé pour chacun des inds
7. si la recherche et fructueuse,
% 1'individus des remplacés par le résultat.

if NSTRAT-=1
for i = l:m 7 boucle sur les individus

7, déterminer le nombre de matériaux différents

'/, et combien de fois il sont représentés

7. vecteurs couplés
matdiff=[]; '/, contient les différents mat du phénotype
combienmatdiff=[]; % contient le nb que fois qu'un mat est présent
unplus=l;
for istrat=l:NSTRAT
lemat=Phenmeil(i,istrat);
yela=find(matdiff==lemat);
if isempty(yela)
matdiff(unplus)=lemat;
combienmatdiff(unplus)=l;
unplus=unplus+l;
else
combienmatdiff(yela)=combienmatdiff(yela)+l;
end
end 7, for istrat

7. mettre dans l'ordre

[matdiff,indice]=sort(matdiff);
UNcombienmatdiff=combienmatdiff;
for iy=l:length(matdiff)
combienmatdiff(iy)=UNcombienmatdiff(indice(iy));
lill.l

'/, c o n s t r u c t i o n des v o i s i n s
y.

'/, récupérer le vecteur porosité

porosite=Phenmeil(i,NSTRAT+l:2*NSTRAT);

"/, nombre de possibilité de combinaison

nbdemat=length(matdiff);
trouvevide=find(matdiff==4) ;

if isempty(trouvevide)
y, c'est ok
elseif NSTRAT-=1
nbdemat=nbdemat-l;
end

if nbdemat==0
display(Phenmeil);
pause(5);
e r r o r C b e n c ' ' e s t vide p a r t o u t ! ! ! ' ) ;
end

'/, i n i t i a l i s a t i o n
un=[];
deux= [] ;
trois=[] ;
Phenlocal=[] ;
Phenlocal=[] ;
Chromlocal= [] ;

'/, vecteur des strates concaténer

if nbdemat==2

un = rep(matdiff(1),[1,combienmatdiff(1)] ) ;
deux = rep(matdiff(2),[1,combienmatdiff(2)]);

'/, combler avec des matériaux vides si trop courts

if "isempty(trouvevide)
manque=NSTRAT-length( [un,deux] ) ; 7, combien de matériaux vide
comblervide=rep(MatDernier, [1,manque] ) ; '/, combler l'espace
PhenlocaKl, :) = [un,deux,comblervide,porosité] ;
Phenlocal(2,:)=[deux,un,comblervide,porosité];
else
P h e n l o c a K l , :) = [ u n , d e u x , p o r o s i t e ] ;
Phenlocal(2,:)=[deux,un,porosité];
end

elseif nbdemat==3
un = rep(matdiff(1),[1,combienmatdiff(1)] );
deux = rep(matdiff(2),[1,combienmatdiff(2)]);
trois = rep(matdiff(3),[1,combienmatdiff(3)]);
 */, combler avec des matériaux vides si trop courts
if "isempty(trouvevide)
7. combien de matériaux vide
manque=NSTRAT-length([un.deux,trois]);
comblervide=rep(MatDernier,[1,manque]);
7, combler l'espace

PhenlocaKl, ) = [un,deux,trois,comblervide, porosité]

Phenlocal(2, ) = [un,trois,deux,comblervide,porosité]
PhenlocalO, )=[deux,un,trois,comblervide,porosité]
Phenlocal(4, )=[deux,trois,un,comblervide,porosité]
Phenlocal(5, )=[trois,un,deux,comblervide.porosité]
Phenlocal(6, )=[trois,deux,un,comblervide,porosité]

else

PhenlocaKl, )=[un,deux,trois,porosité]
Phenlocal(2, ) = [un,trois,deux.porosité]
PhenlocalO, ) = [deux,un,trois.porosité]
Phenlocal(4, )=[deux,trois,un,porosité]
Phenlocal(5, )=[trois,un,deux,porosité]
Phenlocal(6, ) = [trois,deux,un,porosité]
end

end
'/, construction du chromosome
[ Chromlocal ] = rv2bs( Phenlocal, FieldD, ...
PRECI_matID, PRECI_poro, ...
NSTRAT, MatDernier) ;

if nbdemat~=l "/, impossible de faire de combinaisons !

"/. Correction des gènes

[Phenlocal,Chromlocal] = correctgene(Phenlocal,Chromlocal,...
FieldD,PRECI_matID,...
PRECI_poro,NSTRAT,...
MatDernier,porosite_max);

'/, Tronquer le stratifié si l'on rencontre une couche a porosité max.

[Chromlocal,Phenlocal] = tronque(Chromlocal,Phenlocal,...
PRECI_matID,PRECI_poro,...
porosite_max,NSTRAT,MatDernier,...
traduction,Phen_possible);

'/, Evaluation de la population issus de la recherche locale

[QbjVlocal,Dbase,T,v,Prop,strat,algo] = naturelle(strat,algo,...
Phenlocal,NSTRAT,Ra,Dbase,L,contraintes);

'/, Evaluation de la population initiale

[ObjV,Dbase,T,v,Prop,strat,algo] = naturelle(strat,algo,...
PhenmeiKi, :),NSTRAT,Ra,Dbase,L,contraintes) ;

'/, Comparer la population issus de la recherche locale avec l'individus

[Elbest.lequelbest] ■ min(ObjVlocal);
 '/. vérifier si l'amélioration locale est meilleure
if Elbest < ObjV
% l'individu est remplacé
PhenmeiKi, : ) = P h e n l o c a l ( l e q u e l b e s t , :) ;
ChrommeiKi, :)=Chromlocal(lequelbest, : ) ;
ObjV(i)=ObjVlocal(lequelbest);
stumieux=stumieux+l;
V.display(stumieux)
end '/, condition Elbest

end % nbdemat==l
end % boucle en i sur les individus
end */. if NSTRAT ~= 1

algo.stumieux=stumieux;

'/.Fin

B.3.16 Croisement multi-points

'/. xovmp.m - croissement multi-point

y.
'/, Cette fonction prend un matrice de Chromosomes d'individus
°/, en format binaire et procède au croisement des pairs d'individus
'/. consécutif selon une probabilité de croisement Px. Cette fonction
'/, retourne un matrice contenant le résultat des croisements.

"/, Syntaxe: [ NewChrom ] = xovmp(01dChrom, Px, Npt, Rs)

y.
% Paramètre d'entrée:
y.
7. OldChrom Matrice contenant la description des gènes
'/. sous la forme de chromosomes chaque ligne
'/. Correspondant a un individus. Avant le
'/. Croisement.
y. Px Probabilité de croissement
'/. Npt le nombre de points utilisé pour le croisement
y. ( possibilité 1, 2 ou 0-> suffle)
y. Rs Indidique si l'on force la production d'enfant
x Différents de leurs parents. ( 1 ou [])

x
'/, Paramètre de sortie:
y.
'/. NewChrom Matrice contenant la description des gènes
X sous la forme de chromosomes chaque ligne
X Correspondant à un individu. Résultats de du
X Croisement.
X

'X Auteurs: Carlos Fonseca, Andrew Chipperfield

'/, Adaptation: Charles Villemure
7.
y0* *************************************************************************
function NewChrom = xovmpCOldChrom, P x , Npt, R s ) ;
y,**************************************************************************

% identifier la taille de la population

[Nind.Lind] = size(OldChrom);

7. Récupérer les paramètres d'entrée

if Lind < 2, NewChrom = OldChrom; return; end
if nargin < 4, Rs = 0; end
if nargin < 3, Npt = 0 ; Rs = 0; end
if nargin < 2, Px = 0.7; Npt = 0 ; Rs = 0; end
if isnan(Px), Px = 0.7; end '/. probabilité par défaut
if isnan(Npt), Npt = 0; end
if isnan(Rs), Rs = 0; end
if isempty(Px), Px = 0.7; end
if isempty(Npt), Npt = 0; end '/, suffle par défaut
if isempty(Rs), Rs = 0; end "/, non-actif par défaut

Xops = floor(Nind/2);
DoCross = randCXops, 1) < Px; '/, permet déterminer le cross ou pas
odd = l;2:Nind-l;
even = 2:2:Nind;

'/. Calcul la longueur efficace de chaque pair de chromosomes

Mask = "Rs I COldChromCodd, :) "» OldChromCeven, :));
Hask = cumsumCMask')';

'/, Calcul les points de croisement pour chaque pair d'individus, tout
'/, dépendant de leur longueur effective et de la probabilité de croissement
'/, Px C Deux site égaux conduits à aucun croisement)
xsitesC:, 1) = MaskC:, Lind);
if Npt >= 2,
xsitesC:, 1) " ceil(xsites(:, 1) .* randCXops, 1));
end

'/.xsitesC: ,2) = remCxsites + ceil((MaskC : , Lind)-1) .* randCXops, 1)) ...

•/, .* DoCross - 1 , MaskC:, Lind) ) + l;
'/, Adaptation pour éviter des NaN
qtl=xsites + ceilCCMaskC:, Lind)-1) .* randCXops, 1)).* DoCross - 1;
qt2=Mask(:, Lind);
[m,n]=sizeCxsitesC:,1));
for i™l:m
if qtlCi) < 0
qtlCi)=l;
end
if qt2Ci) < 0 | qt2Ci) ==0
qt2Ci)=l;
end
end
xsitesC:,2) = remCqtl,qt2)+l;

'/. Exprime les points de croisement en terme de 0 ou 1 C inactif actif )

Mask = CxsitesC:,ones(l,Lind)) < Mask) == ...
 (xsites(:,2*ones(l,Lind)) < Mask);
if ~Npt,
shuff = rand(Lind,Xops);
[ans.shuff] ■ sort(shuff);
for i=l:Xops
01dChrom(odd(i),:)=01dChrom(odd(i),shuff(:,i));
01dChrom(even(i),:)=01dChrom(even(i),shuff(:,i));
end
end

'/, Croissement
NewChromCodd,:) - (01dChrom(odd,:).* Mask) + (01dChrom(even,:).*(~Mask));
NewChrom(even,:) = (OldChromCodd,:).*(~Mask)) + (OldChromCeven,:).*Mask);

% Si le nombre d'individus est impair, le dernier individu ne peut être

'/. couplé et doit être inclus tel quel dans la nouvelle population,
if remCNind,2),
NewChrom(Nind,:)=01dChrom(Nind,:);
end
if ~Npt,
[ans,unshuff] = sort(shuff);
for i=l:Xops
NewChrom(oddCi),:)=NewChrom(odd(i).unshuffC:,i));
NewChromCevenCi),:)=NewChrom(even(i),unshuffC:,i));
end
end

'/. f i n

B.3.17 Mutation

'/, mut.m - opérateur de mutation sur les gènes

'/.
'/, cette fonction effectue la mutation des gènes des individus d'un
'/, population selon un probabilité donnée.
•/.
7.
'/. Syntaxe: NewChrom = mut(01dChrom,Pm,BaseV)
y.
'/, Paramètre d'entrée:
%
'/. OldChrom - Matrice contenant les chromosomes du
'/, Population. Chaque linge correspond à la chaîne
'/■ de binaires associées à un individu.
'/. C la première variable est à gauche
'/, et la dernière à droite )
y.
'/. Pm - probabilité de mutation
X
'/, BaseV - vecteur contenant la base de chaque variable
'/, C paramètre optionnel, par défaut base =2)
'/, paramètre de sortie:
•/.
7. NewChrom - Matrice conteant la version mutée de OldChrom
•/.
'/.
7,
'/, Auteur: Andrew Chipperfield
'/o Adaptation: Charles Villemure
'/. 20.12.2006
'/.
%***********************++********************+*********+*+****************
function NewChrom ■ mut(OldChrom,Pm,BaseV)
%******************** ***************************************** *************

'/, détermine la dimension de la population

[Nind, Lind] = size(OldChrom) ;

'/, fixe la probabilité de mutation (tenir compte du nombre de bits total)

if nargin < 2, Pm = 0.09 ; end %0.7/Lind
if isnan(Pm), Pm = 0.09 ; end '/,0.7/Lind

"/, Ajuste la base des variables ( par défaut base = 2, binaire)

if (nargin < 3 ) , BaseV = crtbase(Lind); end
if (isnan(BaseV)), BaseV = crtbase(Lind); end
if (isempty(BaseV)), BaseV = crtbase(Lind); end

if (nargin == 3) & (Lind "= length(BaseV))

errorCOldChrom and BaseV are incompatible'), end

% Crée une matrice des marqueurs de mutations

BaseM = BaseV(ones(Nind,1),:) ;

"/, perform mutation on chromosome structure

NewChrom = rem(01dChrom+(rand(Nind,Lind)<Pm)...
.*cell(rand(Nind,Lind).*(BaseM-l)),BaseM);
% fin

B.3.18 Réinsertion

% reins.m _ Réinsertion des enfants dans la population

y.
'/, Cette fonction réinsère les enfants dans la population en remplaçant les
'/. parents.

% Syntaxe: [Chrom, ObjVCh]

"/. = relns(Chrom, SelCh, SUBP0P, InsOpt, ObjVCh, ObjVSel)
y.
'/, Paramètres d'entrée :
y. Chrom - Matrice contenant les chromosomes de la
Z Population. Chaque linge correspond
'/, à la chaine de binaires associés à un individu.
7. (la première variable est à gauche et la dernière à droite)
'/, SelCh - Matrice contenant les enfant s. Chaque linge correspond
'/. à la chaîne de binaires associés à un individu.
'/, (la première variable est à gauche et la dernière à droite)
•/.
•/.
7. Paramètres de sortie :
7. Chrom - Matrice contenant les chromosomes de la nouvelle
% population. Chaque linge correspond
7. à la chaîne de binaires associés à un individu.
7. (la première variable est à gauche et la dernière à droite)
7.
7. Auteur: Hartmut Pohlheim
7e Adaptation: Charles Villemure
7. 15.12.2006

function [Chrom,ObjVCh] = reins(Chrom, SelCh....

SUBPOP, InsOpt, ObjVCh, ObjVSel);
70**************************************************************************

7. récupérer les paramètres

if nargin < 2, error('Paramètres entrés manquants ' ) ; end
if (nargout == 2 & nargin < 6)
errorCDoit contenir: ObjVCh and/or ObjVSel');
end

[NindP, NvarP] = size(Chrom);

[NindO, NvarO] = size(SelCh);

if nargin == 2, SUBPOP - 1; end

if nargin > 2,
if isempty(SUBPOP), SUBPOP = 1;
elseif isnan(SUBPOP), SUBPOP = 1;
elseif length(SUBPOP) "- 1, error('SUBPOP -> un scalaire'); end
end

if (NindP/SUBPOP) '= fix(NindP/SUBPOP)

errorCChrom et SUBPOP ne sont pas de la même taille');
end
if (NindO/SUBPOP) "= fix(NindO/SUBPOP)
errorCSelCh et SUBPOP ne sont pas de la même taille');
end
NIND = NindP/SUBPOP; 7. Calcul le nombre d'individus par sous-pop
NSEL = NindO/SUBPOP; 7. Calcul le nombre d'enfants par sous-pop

IsObjVCh = 0; IsObjVSel = 0;
if nargin > 4,
[mO, nO] = size(ObjVCh);
if nO "■ 1, errorCObjVCh doit être un vecteur colonne'); end
if NindP "- mO
errorCChrom et ObjVCh ne sont pas de la même taille'); 3
end
IsObjVCh = 1;
end
if nargin > 5,
 [mO, nO] = size(ObjVSel);
if nO "■ 1, errorC ObjVSel doit être un vecteur colonne'); end
if NindO "- mO
errorCSelCh et ObjVSel ne sont pas de la même taille');
end
IsObjVSel = 1;
end

if nargin < 4, INSR = 1.0; Select = 0; end

if nargin >= 4,
if isempty(InsDpt), INSR = 1.0; Select = 0;
elseif isnan(InsOpt), INSR = 1.0; Select = 0;
else
INSR = NaN; Select = NaN;
if (length(InsOpt) > 2 ) ,
error('Paramètre InsQpt est trop long');
end
if (length(InsOpt) >= 1 ) , Select = InsOpt(l); end
if (length(InsOpt) >= 2 ) , INSR = Ins0pt(2); end
if isnan(Select), Select = 0; end
if isnan(INSR), INSR =1.0; end
end
end

if (INSR < 0 | INSR > 1)

errorCtaux d insertion doit être un scalaire dans [0, 1]');
end
if (INSR < 1 & IsObjVSel "= 1)
errorC ObjVSel nécessaire pour le choix des enfants');
end
if (Select "= 0 & Select "- 1)
error('Méthode de sélecion doit être 0 ou 1');
end
if (Select == 1 & IsDbjVCh == 0)
errorCObjVCh requis pour les échanges f itness-based ');
end

if INSR == 0, return; end

'/. Nombre d'enfants à insérer
NIns = min(max(floor(INSR*NSEL+.5),l),NIND);

'/, insertion pour chaque sous-pop

for irun = 1:SUBP0P,
7, Calcul de la position dans les aciennes
7, sous-pop où les enfant sont insérés
if Select == 1, % reinsertion fitness-based
[Dummy, Chlx] = sort(-ObjVCh((irun-l)*NIND+l;irun*NIND));
else '/, reinsertion uniforme
[Dummy, Chlx] = sort(rand(NIND,1));
end
Poplx = ChIx((l:NIns)')+ (irun-l)*NIND;
'/, Calculate position of Nins-'/. best offspring
'/. Calcul de la position des meilleurs enfant
if (NIns < NSEL), '/. choix des meilleurs enfants
[Dummy,Offlx] - sort(ObjVSel((irun-l)*NSEL+l:irun*NSEL));
 else
Offlx = (l:NIns)';
end
Sellx = OffIx((l:NIns)')+(irun-l)*NSEL;
'/, insertion des enfant dans la sous-pop — > nouvelle pop
ChromCPopIx,:) = SelChCSelIx,:);
if (IsObjVCh -« 1 k IsObjVSel — 1 ) ,
ObjVCMPopIx) = ObjVSel(Sellx);
end
oilcl

'/. Fin

B.3.19 Sélection des individus

*/, select.m - sélection universelle

*/.
'/, Cette fonction effectue une sélection universelle. La fonction peut
'/. traiter des multi-populations.
y.
% Syntaxe: SelCh = select(SEL_F, Chrom, FitnV, GGAP, SUBPOP)
y.
% Paramètres d'entrée :
'/, SEL_F - Nom de la fonction de sélection
"/, Chrom - Matrice contenant les chromosomes de la
'/. population. Chaque linge correspond
% à la chaîne de binaires associées à un individu.
'/. (la première variable est à gauche et la dernière à droite)
"/, FitnV - Vecteur colonne contenant la valeur associocié à la
"/. performance relative d'un individu.
% GGAP - (optionel) Pourcentage des individus sélectionnés
'/■ SUBPOP - (optionel) nombre de sous populations
'/.
'/. Paramètres de sortie :
'/. SelCh - Matrice contenant les chromosomes de la nouvelle
'/. Population. Chaque linge correspond
'/. à la chaine de binaires associés à un individu.
% (la première variable est à gauche et la dernière à droite)
y.
'/, Author: Hartmut Pohlheim
'/, Adaptation: Charles Villemure
•/. 12.12.2006
%********* + + ************** + ****** + * + + ********* + * * * * * + * ++* + * * + + * * * + * * + * * * + * +
function SelCh = select(SEL_F, Chrom, FitnV, GGAP, SUBPOP);
y**************************************************************************

°/. vérifier les paramètres d'entrée

if nargin < 3, error('paramètres d entrée manquante ' ) ; end

'/, déterminer la taille de la population

[NindCh.Nvar] = size(Chrom);
 [NindF.VarF] = size(FitnV);
if NindCh "- NindF
errorCChrom et FitnV ne sont pas de la même taille');
011(1

if VarF "= 1, errorCFitnV doit être un vecteur colonne ' ) ; end

if nargin < 5, SUBPOP ■ 1; end

if nargin > 4,
if isempty(SUBPOP), SUBPOP = 1;
elseif isnan(SUBPOP), SUBPOP = 1;
elseif length(SUBPOP) ~= 1
errorCSUBPOP doit être un scalaire');
end
end

if (NindCh/SUBPOP) ~= fix(NindCh/SUBPOP)
errorCChrom er SUBPOP ne sont pas de la même taille');
(sud

'/, Déterminer le nombre d'individus par sous-pop

Nind - NindCh/SUBPOP;

if nargin < 4, GGAP = 1; end

if nargin > 3,
if isempty(GGAP), GGAP = 1;
elseif isnan(GGAP), GGAP = 1;
elseif length(GGAP) ~= 1, errorCGGAP soit être un scalaire');
elseif (GGAP < 0 ) , errorCGGAP doit être un scalaire > 0'); end
end

% Calcul le nombre de nouveaux individus à choisir

NSel=max(floor(Nind*GGAP+.5),2);

'/. sélectionne les individus de la pop

SelCh = [] ;
for irun = 1:SUBPOP,
FitnVSub = FitnV((irun-l)*Nind+l:irun*Nind);
ChrIx=feval(SEL_F, FitnVSub, NSel)+(irun-l)*Nind;
SelCh=[SelCh; Chrom(ChrIx,:)];
end

*/. fin