Identification des Paramètres Mécaniques
Identification des Paramètres Mécaniques
net/publication/30517879
Article
Source: OAI
CITATIONS READS
9 202
1 author:
Kossi Atchonouglo
University of Lomé
50 PUBLICATIONS 218 CITATIONS
SEE PROFILE
All content following this page was uploaded by Kossi Atchonouglo on 19 December 2013.
THESE
pour l’obtension du Grade de
DOCTEUR DE L’UNIVERSITE DE POITIERS
(Faculté des Sciences Fondamentales et Appliquées)
(Diplôme National - Arrêté du 7 août 2006)
Ecole Doctorale Sciences pour l’Ingénieur & Aéronautique
Spécialité : Mécanique des Solides, des Matériaux, des Structures et des Surfaces
Claude VALLEE
Directeurs de Thèse :
Magolmèèna BANNA
JURY
Cette thèse a été effectuée en grande partie au Laboratoire de Mécanique des Solides de l’Uni-
versité de Poitiers, dirigé par Monsieur le Professeur Olivier Bonneau. Je lui suis très reconnaissant
pour tous ses appuis qui ont facilités mes conditions de travail.
Je remercie vivement tout particulièrement Monsieur le Professeur Claude Vallée, qui m’a
proposé ce sujet de recherche. Je lui exprime ma gratitude pour sa confiance, sa permanente
disponibilité et son soutien qui m’ont été d’un encouragement et une aide pertinente tout au long
de mes travaux. Il est, à mon avis, le directeur de thèse rêvé de tout doctorant !
Mes remerciements vont également à Monsieur le Professeur Fabrice Bremand, qui m’a accueilli
au sein de l’Equipe Photomécanique et Rhéologie et m’a permis de participer à plusieurs congrès.
Très sensible à l’attention qu’ils ont déjà apporté à mes travaux et à l’honneur que me font
Messieurs Zhi-Qhiang Feng, Professeur à l’Université d’Evry et Kossi Napo, Professeur à l’Univer-
sité de Lomé, en acceptant de juger mes travaux et d’en être rapporteurs, je leur exprime ma très
respectueuse gratitude. Le Professeur Napo s’est déjà illustré à plusieurs reprises par un soutien
sans faille et une sollicitude à toute épreuve lors de mes séjours à l’Université de Lomé.
Que Madame Danielle Fortuné d’une part, Maître de Conférences à l’université de Poitiers, qui
a été ma tutrice pédagogique en enseignement et qui m’a guidé pendant mes travaux et Monsieur le
Professeur Valérie Valle, qui a pu me faire passer de théoricien en praticien d’autre part, trouvent
vi
Une mention spéciale pour Monsieur le Professeur Saïd Zeghloul, aux informaticiens Franck et
Mathieu, Jean-Jacques de la reprographie, Madame Françoise Foures de la Bibliothèque SP2MI,
Djam, Lolotte et Nath de l’administration ; ils savent à quel point le présent travail doit à leurs
efforts.
Mes remerciements vont à tous les membres du Laboratoire de Mécanique des Solides de
l’Université de Poitiers et du Laboratoire sur l’Energie Solaire de l’Université de Lomé pour la
contribution que chacun d’eux m’a apporté.
L’esprit d’équipe entre les doctorants du laboratoire a favorisé considérablement des échanges
très constructifs. Sincères remerciements à Arnaud, Armel, Camélia, Claude, Eric, Khalil, Inès,
Jean-Claude, Lilia, Marie-Eve, Mathieu, Paul, Stephen, Tony, ...
Mes remerciements vont enfin à ma famille et à toutes ces personnes chères à mon cœur (Atchez,
Christophe, Verdier,...) qui n’ont ménagé aucun effort pour m’apporter leur soutien moral constant,
leur appui logistique et financier, sans lesquels ce travail n’aurait pas été possible. Puisse Dieu les
combler de sa grâce et rendre à chacun selon son mérite.
TABLE DES MATIÈRES vii
Notations et abréviations xv
Introduction Générale 1
Introduction Générale 1
Problématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Plan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Bibliographie 146
II.1 Poutre étirée par une force P et soumise à une force transversale . . . 34
II.2 Génération de fonctions linéaires par morceaux . . . . . . . . . . . . . 41
t : temps
x : coordonné cartésienne
xi : coordonnés cartésiennes (i = 1, 2, 3)
y : coordonné cartésienne
z : coordonné cartésienne
φi : Fonction élémentaires
ζ : pas de discrétisation spatiale
Mécanique
a(M ) : Vecteur accélération du solide
F : Forces extérieures à l’exclusion du poids du solide
I : Matrice inertielle classique
J 0 : Matrice d’inertie de Poinsot généralisée calculée à l’instant initial
J : Matrice d’inertie de Poinsot généralisé
Mex : Moment des forces extérieures
m : masse
R : Matrice rotation 3 × 3
R0 : Repère fixe de R3
R00 : Repère fixe de R4
R : Repère lié au solide de R3
R0 : Repère lié au solide de R4
T : Vecteur translation
V(M ) : Vecteur vitesse d’un point M du solide
xvi
χ : Résultante dynamique
Υ : Torseur des forces et des moments : matrice 4 × 4
contenant les composantes des forces extérieures et des moments
Θ : Matrice des accélérations dans R00
Ω : Matrice des rotations 4 × 4
δ : Moment dynamique
Thermique
c : Chaleur spécifique
C : Matrice de Capacité
e : Energie interne spécifique
E : Energie interne spécifique
F : Vecteur résultant des conditions aux limites
K : Matrice de conductivité
h : Coefficient d’échange de chaleur par convection
hr : Coefficient d’échange de chaleur par rayonnement
k : Conductivité thermique
−
→
n : Vecteur tnormal unitaire, orienté vers l’extérieur
q : Flux de chaleur constant
s : Entropie sécifique
S : Entropie
T (x, t) : fonction température
Ta : Température ambiante
TΣ : Température de la surface Σ
(α0 , ..., αn ) : n + 1 variables d’état définissant les propriétés du système
Ω : Domaine géométrique occupé par le matériau
Σ : Bord de Ω
σ : Constante de Stefan-Boltzmann
: Emissivité : facteur de forme gris
ψ : Flux thermique
Notation xvii
Problématique
Les lois physiques reliant la température aux flux de chaleur ou les déformations
aux contraintes font intervenir des paramètres thermiques et mécaniques (conducti-
vité thermique, module d’Young, etc,...) qui dépendent de la nature des matériaux ;
ce sont les lois de comportement. Si les gradients de température ou les déforma-
tions sont petites, les lois de comportement sont en général linéaires. De ces lois de
comportement et des bilans énergétiques, on déduit des équations d’évolution. Ces
équations d’évolution sont décrites par des équations différentielles ordinaires (EDO)
ou aux dérivées partielles (EDP). La linéarité des lois de comportement conduit à
la linéarité des équations d’évolution. Les équations aux dérivées partielles linéaires
que nous aborderons sont du type
∂2y ∂y ∂2y ∂y
a + b + c + d = f (x, t). (1)
∂x2 ∂x ∂t2 ∂t
Les constantes a, b, c et d caractérisent le matériau. Les fonctions y et f dépendent
des variables d’espace x et du temps t.
Chercher à exprimer y en fonction de x et t connaissant les constantes a, b, c, d
et la fonction f c’est résoudre le problème direct.
Le problème inverse se pose de deux manières.
Une première branche porte sur la recherche de la fonction f qui conduit pour des
constantes a, b, c et d connues, à une fonction y souhaitée. C’est le cas de la recherche
du flux de chaleur conduisant à un champ de température souhaité.
Une seconde branche concerne la recherche des constantes a, b, c et d lorsque y et f
sont connues. Par exemple, on cherche la valeur de la conductivité thermique d’un
matériau, le flux de chaleur et les températures étant connues.
Le but de notre travail est de développer des approches numériques permettant
de résoudre les problèmes de la seconde branche. Ces approches sont appliquées à
2 Introduction
une situation mécanique et à une situation thermique faisant intervenir des équa-
tions du type (1).
La situation mécanique concerne l’identification des dix caractéristiques d’inertie
(masse, position du centre d’inertie et tenseur d’inertie) d’un solide rigide.
Dans l’analyse du mouvement d’un solide rigide, le calcul de la quantité de mouve-
ment et du moment cinétique nécessite l’identification préalable de ses caractéris-
tiques inertielles. De plus, en statique, un système est en équilibre si la projection
de son centre d’inertie sur le sol est dans son polygone de sustentation. La finalité
de ces recherches est l’identification des caractéristiques d’inertie d’un athlète dans
le souci d’améliorer ses performances. En supposant l’athlète composé de 16 troncs
de cône rigides, il s’agit donc d’identifier 160 caractéristiques d’inertie. C’est dans
ce cadre que j’ai effectué des travaux de recherche en collaboration avec l’équipe du
Geste Sportif de notre Laboratoire.
La situation thermique concerne les études thermomécaniques des matériaux,
spécialement les polymères. Ces études, réalisées au sein de l’équipe Photoméca-
nique et Rhéologie dont je fais partie, nécessitent la connaissance préalable de la
conductivité thermique et de la chaleur volumique des matériaux. Notre travail est
de développer et de valider une nouvelle méthode d’identification des paramètres
thermophysiques. Pour cela, nous avons élaboré une méthode d’optimisation globale
dans le cas unidimensionnel. Les paramètres thermophysiques sont identifiés à partir
des champs de température obtenus expérimentalement pour des flux de chaleur de
chauffe fixés. Ce travail est lié aux recherches menées au « Laboratoire sur l’Ener-
gie Solaire de l’Université de Lomé » sur la caractérisation de capteurs hybrides
photovoltaïques intégrés aux bâtiments.
Une étape préliminaire à la résolution du problème inverse est la modélisation du
phénomène : c’est le problème direct qui décrit comment les paramètres du modèle
se traduisent en effets observables expérimentalement. Ensuite, à partir de mesures
expérimentales, la démarche va consister à approcher au mieux les paramètres qui
permettent de rendre compte de ces mesures.
Avant d’envisager la résolution d’un système d’équations aux dérivées partielles
c’est-à-dire la résolution du problème direct, il convient d’analyser les équations
du modèle. À ce niveau, on doit se poser des questions sur l’existence, l’unicité
de la solution, la sensibilité des solutions aux perturbations, la croissance ou la
décroissance des solutions en fonction du temps, etc. Cette analyse conduit à la
résolution de problèmes mathématiques extrêmement complexes, qui cependant sont
Introduction 3
Plan
Les algorithmes élaborés au chapitre I sont validés dans la suite en les appliquant à
deux situations mécaniques :
1. chapitre II : Application à l’identification de la matrice tridiagonale symétrique
issue de la modélisation par élements finis d’une équation différentielles linéaire
du second ordre à coefficients variables. Il s’agit d’identifier les coefficients de
cette équation différentielle. Les coefficients variables se retrouvent incorporés
dans des intégrales permettant de calculer les coefficients d’une matrice qui
se révèle symétrique tridiagonale. Dans une première étape, les valeurs appro-
chées des coefficients de la matrice symétrique tridiagonale sont identifiées par
application de l’algorithme élaborés au chapitre I. Nous développons ensuite un
algorithme permettant d’en déduire une approximation linéaire par morceaux
des coefficients variables de l’équation différentielle de départ.
2. chapitre III : Identification des dix caractéristiques d’inertie d’un solide rigide.
Nous regroupons les composantes du moment dynamique et de la résultante
dynamique dans une matrice 4 × 4 antisymétrique. Nous faisons de même pour
la résultante dynamique et le moment résultant des efforts extérieurs. Les
caractéristiques d’inertie du solide se retrouvent dans une matrice 4 × 4 symé-
trique défini-positive J . L’équation du mouvement est ramenée à une équation
dans l’espace des matrices 4 × 4 antisymétriques. Elle exprime l’égalité entre
le torseur dynamique et le torseur des efforts extérieurs sous une forme ma-
tricielle appropriée à l’identification de la matrice J . Une formulation directe
des équations du mouvement du solide rigide dans l’espace des matrices 4 × 4
est établie par application des travaux virtuels. L’identification de la matrice
J restitue la masse, la position du centre d’inertie et le tenseur d’inertie du
solide.
Dans la deuxième partie, nous commençons par rappeler au chapitre IV la fromu-
lation de l’équation aux dérivées partielles de Fourier régissant la conduction de la
chaleur dans un solide. La discrétisation spatiale par éléments finis de l’équation aux
dérivées partielles nous conduit à remplacer le problème direct par la résolution d’une
équation différentielle linéaire. C’est cette équation qui sera utilisée pour résoudre
les problèmes inverses, c’est-à-dire pour identifier paramètres thermophysiques mis
en jeu (conductivité thermique k et chaleur volumique ρc).
Au chapitre V, nous avons développé une nouvelle technique d’identification si-
multanée de la conductivité thermique k et de la chaleur volumique ρc. Une première
6 Introduction
Première partie
Sommaire
I.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 10
I.2 Position du problème . . . . . . . . . . . . . . . . . . . . . 11
I.3 Définitions et notations en algèbre linéaire . . . . . . . . 11
I.4 Identification de la matrice A . . . . . . . . . . . . . . . . 14
I.5 Identification de A symétrique . . . . . . . . . . . . . . . . 17
I.5.1 Multiplicateurs de Lagrange . . . . . . . . . . . . . . . . 17
I.5.2 Elaboration de l’algorithme d’identification de A symétrique 19
I.5.3 Algorithme de projection pour les matrices symétriques . 21
I.6 Identification de A symétrique tridiagonale . . . . . . . . 21
I.6.1 Elaboration de la procédure d’identification de A symé-
trique tridiagonale . . . . . . . . . . . . . . . . . . . . . . 24
I.6.2 Algorithme de la Méthode de Projection des Matrices sy-
métriques Tridiagonales . . . . . . . . . . . . . . . . . . . 26
I.7 Méthode intermédiaire . . . . . . . . . . . . . . . . . . . . 28
I.7.1 Elaboration de la méthode d’identification . . . . . . . . . 28
I.7.2 Algorithme d’identification . . . . . . . . . . . . . . . . . 31
I.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
10 I. Un problème Type d’Identification
I.1 Introduction
D’une manière générale, les matrices issues de la discrétisation par éléments finis
ou différences finies mettent en jeu les caractéristiques mécaniques du matériau.
Les matrices symétriques tridiagonales sont aussi très courantes dans la résolution
numériques des problèmes du type Sturm-Liouville [17].
Rappelons que l’espace des matrices symétriques et l’espace des matrices tridia-
gonales et symétriques sont deux sous-espaces vectoriels de l’espace vectoriel des
matrices carrées.
Dans le cas du système (I.1), le problème direct consiste à déterminer le vecteur X
pour A et b connus. Dans la pratique, il s’agit de déterminer le moment fléchissant
d’une poutre connaissant les caractéristiques d’inertie et les forces. Le problème
inverse consiste à identifier la matrice A de l’équation (I.1). Lorsque cette matrice
est de faible taille, les méthodes habituelles de résolution de système linéaire sont
bien adaptées. Pour des matrices de grande taille ( matrices de taille supérieure à
4), certaines méthodes restituent difficilement la matrice attendue.
Le but de ce chapitre est de tenir compte de la symétrie et de la tridiagonalité
de la matrice A lors de la résolution numérique et de chercher à réduire le coût des
opérations.
Nous décrirons des procédure d’identification basée sur la la méthode des moindres
carrés tenant compte de la symétrie ou de la tridiagonalité de A par introduction
I.2. POSITION DU PROBLÈME 11
Proposition I.3.1
Pour u et v deux vecteurs de Rn , on a :
u ⊗ v = uvT
Proposition I.3.2
Soient u et v deux vecteurs de Rn . Le produit scalaire de u et v est le réel
hu, vi = tr(u ⊗ v)
Proposition I.3.3
Soient A une matrice de Mn et deux vecteurs u et v de Rn . Le produit scalaire de
Au et v est
hAu, vi = tr(A(u ⊗ v)).
Proposition I.3.4
Soient P et Q deux matrices n × n :
(i) le produit scalaire est symétrique : tr(P Q)= tr(QP ) ;
(ii) une matrice et sa transposée ont la même trace : tr(P T )= tr(P )
et
n
X
UV = λi (ui ⊗ ui )V
i=1
Xn
= λi ui ⊗ (V ui )
i=1
Par suite
n
X
tr(U V ) = λi hui , V ui i.
i=1
Nous allons élaborer une méthode pour identifier la matrice A du système li-
néaire :
AX = b (I.3)
lorsque les vecteurs X et b sont connus. Plus précisement, nous allons construire
une procédure pour résoudre le problème inverse suivant :
"plusieurs couples (X(i) , b(i) ), indescés de i = 1 à m ont été mesurés, peut-on iden-
tifier l’application linéaire A ?
Notre procédure d’identification est basée sur la méthode des moindres carrés
Pour déterminer la matrice A, nous minimiserons donc la forme quadratique
m
1X
F (A) = kAX(i) − b(i) k2 (I.4)
2 i=1
Les couples (X(i) , b(i) ) satisfont AX(i) = b(i) et m est le nombre de couples connus.
La matrice A sera identifiée par annulation du gradient de la forme quadratique
F , i.e. en résolvant le système linéaire
−−→
gradF (A) = 0
Il nous reste à calculer la dérivée de la forme quadratique F (A). Pour cela, il est utile
de rendre euclidien l’espace des matrices n × n en le munissant du produit scalaire
des matrices (I.3.7). La généralisation du gradient est alors
∂F
[(gradF )(A)]ij =
∂Aij
Remarque I.4.1
Pour calculer la dérivée de la forme quadratique F , nous allons utiliser la remarque
suivante de Jean-Marie Souriau.
F (A + αδA) − F (A)
lim = hgradF (A), δAi
α→0 α
Annuler le gradient de F en A revient à annuler la dérivée par rapport à α de
F (A + αA) pour toute variation δA.
16 I. Un problème Type d’Identification
+ 21 tr AT AH − AB T + 12 m
P (i) 2
i=1 kb k
Les matrices H et δA(δA)T sont symétriques et positives. Le réel tr H(δA)T δA
est donc positif d’après la proposition (I.3.5). Par suite le trinôme F (A + αδA) est
convexe en α.
La dérivée de F (A + αδA) en α égale à zéro vaut
AH = B (I.6)
Dans l’espace des matrices carrées, le gradient de la forme quadratique F (A) est[96]
∂F
= (AH − B)ij .
∂Aij
Remarque I.4.2
Nous pouvons obtenir directement l’équation (I.6) en multipliant AX(i) = b(i) ten-
soriellement par le vecteur X(i) puis en sommant sur le nombre de mesures :
m
X m
X
(i) (i)
A X ⊗X = b(i) ⊗ X(i) (I.7)
i=1 i=1
AH = B
Remarque I.4.3
Bien qu’elle soit basée sur la méthode des moindres carrés, notre procédure d’identi-
fication se révèlera donc exacte lorsque la matrice H sera inversible.
Remarque I.4.4
Si m est supérieur ou égal à n, la matrice H est attendue défini-positive, et donc
inversible [99]. Et expérimentalement ceci revient à faire plusieurs mesures "indé-
pendantes". Pour une poutre en flexion pure par exemple, il faut plusieurs mesures
du moment fléchissant pour des forces distinctes.
Nous allons forcer la matrice A à être symétrique. Il suffit pour cela de reprendre
la minimisation de (I.4) en respectant les conditions linéaires
A = AT
−−→
Lagrange (ΛT = −Λ). La variation DF (A)δA = hgradF (A), δAi, δA ∈ Mn ne doit
plus être nulle mais égale à hδf (A), Λi.
1 T
hδf (A), Λ = tr Λ (δA − δA )
2
1 1
= tr(ΛδA) − tr(ΛδAT )
2 2
= tr(δAΛ)
Ainsi on conclura à
AH = B + Λ (I.8)
Commentaire I.5.1
La méthode que nous développons ici est est inspirée de la technique de
minimisation avec des contraintes, à la différence que nous restons dans
l’espace des matrices.
On pourrait chercher à déterminer la matrice symétrique A directe-
ment par la minimisation de la fonction
m n
1X X
F (A) = kAXi − bi k + λj (aij − aji )2 ,
2 i=1 j=1
soit n(n + 1)/2 équations permettant de déterminer les n(n + 1)/2 composantes
indépendantes de la matrice symétrique A.
Sur la base des matrices symétriques, A peut se décomposer :
n(n+1)/2
X
A= ak S k ,
k=1
ak , k = 1, ..., n(n + 1)/2 étant des réels. En fait ce sont les composantes de la matrice
A.
Remplaçons dans l’équation de projection, la matrice A par sa décomposition
sur la base Sj ; alors l’identification de la matrice A revient à la détermination des
réels ak dans
n(n+1)/2
X
ak tr(Sj Sk H) = tr(Sj B) j = 1, ..., n(n + 1)/2. (I.9)
k=1
Mjk = tr(Sj Sk H)
Remarque I.5.1
Toute matrice symétrique et positive de Mn est toujours diagonalisable et admet n
valeurs propres réelles positives.
P 2
Cette expression de j a j S j est injectée dans celle de hv, M vi :
XX n
X
aj ak tr(Sj Sk H) = λl tr(Vl VlT H)
j k l=1
Xn
= hλl Vl , HVl i
l=1
C’est positif en tant que somme de nombres positifs puisque H est positive.
Puisque H est défini-positive, la somme ne peut être nulle que si tous les λl sont
nuls. Et dans ce cas puisque H est défini-positive on aurait ( j aj Sj )2 = 0, c’est-
P
P
à-dire j aj Sj = 0. Et comme les Sj constituent la base de l’espace des matrices
symétriques, donc les aj sont nuls.
La matrice du système linéaire de composantes tr(Sj Sk H) est défini-positive.
La résolution de (I.9) tient bien compte de la symétrie imposée sur A.
Ainsi donc pour identifier la matrice symétrique A connaissant plusieurs couples
(X, b), il revient à résoudre le système linéaire AH = B + Λ, où Λ est une ma-
trice antisymétrique. La résolution se fait par projection sur l’espace des matrices
symétriques.
I.6. IDENTIFICATION DE A SYMÉTRIQUE TRIDIAGONALE 21
zj = tr(Sj B) ;
Remarque I.5.2
On vient de résoudre l’équation AH = B+Λ pour A inconnue. La différence est qu’en
utilisant cet algorithme solution, on est sûr d’avoir la matrice A dans le sous-espace
des matrices symétriques.
−β1 λ1 × ··· · ··· ×
× β1 λ1 − β2 λ2 × · ··· ··· ×
. . . . . ..
.. .. .. .. ..
× .
. ..
. . . . .
.. .. .. .. .. .. .
.
.. . .. . .. . .. . .. . .. ..
tr(AΛ) = tr
.
. . . . . . ..
.. .. .. .. .. ..
.
.
.. . .. . .. . .. . .. . ..
×
× ··· · · × × −βn λn + βn−1 λn−1
× ··· ··· ··· × × βn λn
Signalons que nous n’avons représenté que les termes diagonaux du produit AΛ
indispensables pour le calcul de la trace.
c’est-à-dire
β = (Aii+1 + Ai+1i )/2
i
λi = (Aii+1 − Ai+1i )/2.
La décomposition est orthogonale car tr(ΛA) = 0.
c’est-à-dire
tr(δA(HAT − B T ) = tr(δAΛT )
ou encore
AH = B + Λ. (I.12)
I.6. IDENTIFICATION DE A SYMÉTRIQUE TRIDIAGONALE 25
où les γp sont des réels. Ce sont les composantes de la matrice symétrique tridiagonale
A. multiplions l’équation (I.12) par Tm et prennons la trace :
Proposition I.6.2
La matrice symétrique M de composantes Mmp = tr(Tm Tp H) du système (I.14) est
défini-positive.
Preuve. En effet
2n−1
" ! ! #
X 2n−1
X X X
γp γp tr(Tm Tp H) = tr γm Tm γp Tp H
m=1 p=1 m p
!2
X
= tr γm Tm H
m
La matrice symétrique !2
X
γm Tm
m
26 I. Un problème Type d’Identification
Ainsi
XX n
X
γm γp tr(Tm Tp H) = µl tr(Wl WlT H)
m p l=1
Xn
= µl (WlT HWl )
l=1
zj = tr(Tj B)
Exemple I.6.1
Des couples (Xi , bi ), i = 1, ..., 4 sont obtenus d’un système linéaire AX = b, où A
est une matrice 4 × 4 symétrique tridiagonale.
Notons par
−20, 3328 −16, 6024 0, 630162 −25, 4173
−31, 2656 −24, 5509 −1, 53237 −16, 4179
c = X1 X2 X3 X4 =
−31, 7648 −19, 4545 −3, 58713 0, 456126
−21, 89 −11, 054 −0, 927636 −0, 590235
la matrice obtenue par les vecteurs X
200 200 −100 1100
400
400 40 0
d = b1 b2 b3 b4 =
600 400 200 −400
800 200 −100 100
la matrice obtenue par les vecteurs b.
AH = B + Λ
A = (B + Λ)H −1 .
I.7. MÉTHODE INTERMÉDIAIRE 29
AT = A =⇒ H −1 (B T − Λ) = (B + Λ)H −1
=⇒ (B T − Λ)H = H(B + Λ) en multipliant à gauche et à droite par H
HΛ + ΛH = B T H − HB. (I.16)
AH = B + Λ
avec Λ solution de
HΛ + ΛH = B T H − HB.
En procédant de cette manière, la symétrie de A est assuré.
En effet, d’une part, nous avons
Et comme AT = A, donc
AH + HAT = B + B T (1)
on obtient
LH = HΛ + ΛH.
L’équation (I.16) est résoluble si l’application LH est inversible. Ceci peut se vérifier
par le calcul du déterminent de LH .
En effet, soit {ei , i = 1, ..., n} une base de vecteurs propres de la matrice H
choisie orthonormée (l’existence de cette base est assurée par la défini-positivité de
la matrice H).
Choisissons ensuite dans E la base orthonormée (i < j par exemple)
1
√ (ei ⊗ ej − ej ⊗ ei )
2
Dans cette base nous déterminons la matrice de LH par le calcul des images des
éléments de la base de E :
1 1
LH √ (ei ⊗ ej − ej ⊗ ei ) = √ [(Hei ) ⊗ ej − (Hej ) ⊗ ei + ei ⊗ (Hei ) − ej ⊗ (Hei )]
2 2
1
= √ (hj + hi )(ei ⊗ ej − ej ⊗ ei )
2
où hi (resp. hj ) est la valeur propre de la matrice H correspondant au vecteur propre
ei (resp. ej ).
Alors dans la base choisie, la matrice de Lh est la matrice diagonale de composantes
diagonales hi + hj .
Montrons à présent que l’application LH est auto-adjointe.
Soit W une autre matrice antisymétrique de E ; l’espace E est rendu euclidien par
le produit scalaire
hΛ, W i = −tr(HW ).
Calculons alors hLH Λ, W i :
Ceci prouve que l’application LH est auto-adjointe dans E. C’est pourquoi la mé-
thode du gradient conjugué peut s’appliquer. La matrice de LH est diagonalisable
puisqu’elle est auto-adjointe ; il existe donc une base où LH est diagonalisable.
Les valeurs propres de la matrice de LH sont les (n−2)(n−1)/2 réels positifs hi +hj .
Les réels hi + hj sont strictement positifs si la matrice H est défini-positive.
Q
le réel détLH = i<j (hi + hj ) étant strictement positif, l’application LH est donc
auto-adjointe et positive dans l’espace E.
On peut donc déterminer la matrice Λ de l’équation (I.16).
HΛ + ΛH = B T H − HB
AH = B + Λ
I.8 Conclusion
Sommaire
II.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 33
II.2 Problème direct . . . . . . . . . . . . . . . . . . . . . . . . 35
II.2.1 Formulation variationnelle . . . . . . . . . . . . . . . . . . 36
II.2.2 Dérivée de Gateaux de la fonctionnelle J(v) . . . . . . . . 38
II.2.3 Formulation faible . . . . . . . . . . . . . . . . . . . . . . 39
II.2.4 Résolution du problème direct . . . . . . . . . . . . . . . . 41
II.3 Problème inverse . . . . . . . . . . . . . . . . . . . . . . . . 44
II.3.1 Hypothèse de détermination de la matrice A . . . . . . . . 44
II.3.2 Procédure de détermination de la matrice A . . . . . . . . 44
II.3.3 Procédure d’identification de p et de q . . . . . . . . . . . 45
II.4 Application à un cas simple . . . . . . . . . . . . . . . . . 48
II.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
II.1 Introduction
Une équation différentielle est une relation entre une ou plusieurs fonctions incon-
nues et leurs dérivées. Les équations différentielles interviennent dans la modélisation
34 II. Application à l’Identification des Coefficients d’une Equation Différentielle Linéaire
Fig. II.1 – Poutre étirée par une force P et soumise à une force transversale
L’espace Cc∞ (Ω) est l’espace des fonctions de classe C ∞ à support compact dans Ω.
Muni de l’application
Z 21
1 0 2 2
v ∈ H (Ω) −→ kvkH 1 = |v | + |v| )dx ,
Ω
Définition II.2.4 Soit Cc∞ (Ω) l’espace des fonctions de classe C ∞ à support com-
pact dans Ω. L’espace de Sobolev H01 (Ω) est défini comme l’adhérence de Cc∞ (Ω)
dans H 1 (Ω).
Définition II.2.5 Pour un entier m ≥ 0, l’espace de Sobolev H m (Ω) est défini par
où la dérivée partielle
∂αv ∂ α1 +...αn
∂αv = = v, α1 + ...αn = α
∂xα ∂xα1 1 ...∂xαnn
où Π est une fonction continue du triplet (x, v, v 0 ) ainsi que ses dérivées première et
seconde par rapport à x, v et.
Supposons que la fonction v soit continue, qu’elle possède une dérivée continue
0
v sur ]0, L[ et qu’elle prenne aux bornes de l’intervalle [0, L] les valeurs données
Dans ces conditions, une fonction u minimisant J(v) satisafait l’équation d’Euler-
Lagrange
d ∂Π ∂Π
− = 0. (II.8)
dx ∂v 0 ∂v
Nous supposons que u satisfait les conditions nécessaires (de régularité) pour réaliser
le minimum de J(v) [17] et [69].
Le théorème d’Euler-Lagrange assure que si u minimise J(v), alors elle est solu-
tion de l’équation (II.3).
Théorème II.2.1
Si u minimise la fonctionnelle J(v) alors u est solution de l’équation différentielle
(II.3)
Il est facile de voir que le système (II.8) coïncide avec le système (II.3) en choisissant
1 2 1
Π(x, v(x), v 0 (x)) = (p(x)v 0 (x)) + q(x) (v(x))2 − f (x)v(x).
2 2
Puisque la fonction v est supposée continue et dérivable de dérivée v 0 continue, alors
le problème direct peut donc être formulé de manière équivalente sous la forme
variationnelle :
Minimiser dans H 1 ([0, L]) sous la condition v(0) = 0 et v(L) = 0 la fonctionnelle
[17], [29] et [30]
Z L
1 0 2 1 2
J(v) = (p(x)v (x)) + q(x) (v(x)) − f (x)v(x) dx (II.9)
0 2 2
les fonctions p, q et f sont respectivement dans C 1 ([0, L]), C([0, L]) et L2 ([0, L]).
Elles sont données ; la fonction p est strictement positive et la fonction q est positive
et nulle à l’origine ; la fonction p peut être regardée comme un élément de H 1 ([0, L]).
38 II. Application à l’Identification des Coefficients d’une Equation Différentielle Linéaire
v ∈ V −→ J 0 (u, v)
(II.3) par v ∈ H01 ([0, L]) suivie d’une intégration par partie :
Z L
((p(x)u0 )0 (x)v(x) + q(x)u(x)v(x) − f (x)v(x))dx
0
Z L
= (p(x)u0 (x)v 0 (x) + q(x)u(x)v(x) − f (x)v(x))dx − (v(L)u0 (L) − v(0)u0 (0))
0
Z L
= (p(x)u0 (x)v 0 (x) + q(x)u(x)v(x) − f (x)v(x))dx car v ∈ V = H01 ([0, L])
0
Il peut y arriver qu’il n’y ait pas de formulation variationnelle du problème direct,
et dans ce cas, on procède par la formulation faible.
Nous résolvons le problème direct avec les conditions u(0) = u(L) = 0 ; pour
se ramener aux conditions u(0) = ξ, u(L) = η, il suffit de faire un changement de
fonction inconnue en soustrayant la fonction ξ(L − x) + ηx.
On construit Vζ de manière à ce qu’il soit une bonne approximation de H 1 ([0, L]) et
que la solution uζ dans Vζ de la formulation faible soit la plus "proche" possible de
la solution exacte u dans H 1 ([0, L]).
Pour ce faire, on découpe l’intervalle [0, L] en n morceaux égaux de longueur ζ = Ln
(n est supposé supérieur ou égal à 1).
Le sous-espace Vζ est engendré par les fonctions de base φi , i = 1, ..., n [58] et
[87] :
1
ζ x + 1 − i si (i − 1)ζ ≤ x ≤ iζ
φi (x) = − ζ1 x + 1 + i si iζ ≤ x ≤ (i + 1)ζ
0
sinon.
Les fonctions φi sont dans L2 ([0, L]), leurs dérivées au sens des distributions sont
dans L2 ([0, L]), donc les fonctions φi sont dans H 1 ([0, L]).
Le sous-espace Vζ de base{φi , i = 1, ..., n} est donc un sous-espace de dimension
finie de H 1 ([0, L]), sa dimension est n.
Soit uζ ∈ Vζ , alors elle est définie de manière unique par ses valeurs cj = uζ (xj ) aux
n nœuds xj = jζ, 1 ≤ j ≤ n
n
X
uζ (x) = cj φj (x). (II.15)
i=1
42 II. Application à l’Identification des Coefficients d’une Equation Différentielle Linéaire
La méthode des éléments finis consiste à remplacer l’équation (II.14) par sa projec-
tion sur le sous-espace Vh [30]
Théorème II.2.3
On suppose la fonction q positive et p strictement positive.
(1) Étant donné un sous-espace Vζ de dimension finie de l’espace H 1 ([0, L]), il existe
un et un seul élément uζ ∈ Vζ qui vérifie
(2) Cet élément est également caractérisé comme la solution unique du problème :
trouver uζ ∈ Vζ tel que
J(uζ ) = inf J(vζ ).
vζ ∈Vζ
Proposition II.2.1
La matrice n × n A de composantes
Z L
Aij = (Aφi , φj ) = [p(x)φ0i (x)φ0j + q(x)φi (x)φj (x)]dx (II.18)
0
Ce n’est possible que si tous les coefficients ci sont nuls. La matrice A est donc
défini-positive.
Commentaire II.2.1
La matrice A hérite sur le sous-espace Vζ de la symétrie et de la positivité
de A. La coercivité de A entraîne la défini-positivité de A.
Posons Z L
bi = f (x)φi (x)dx (II.19)
0
et interprétons les ci comme valeurs de la fonction uζ aux nœuds d’abscisse xi .
Les valeurs optimales des composantes ci sont obtenues en résolvant le système
linéaire X
Aij cj = bi .
j
Ac = b. (II.20)
44 II. Application à l’Identification des Coefficients d’une Equation Différentielle Linéaire
Connaissant les couples (u(x), f (x)), peut-on trouver deux fonctions p stricte-
ment positive et q positive telles que l’équation différentielle (II.3) soit vérifiée.
La procédure que nous allons suivre consiste à déterminer l’application linéaire
A ou sa version discrétisée A à partir de plusieurs données de couples (uk , F k )
solutions de (II.16). Ensuite les informations sur les fonctions p et q seront déduites
de la connaissance des composantes de la matrice A.
Connaissant les solutions uk (x) pour diverses excitations f k (x), nous souhaitons
en extraire des informations sur les fonctions p et q.
Dans une situation expérimentale, les f k (x) peuvent être les forces exercées sur
une poutre, les uk (x) sont les flèches. Au cours de m expériences indexées par k,
f k (x) et uk (x) sont mesurées aux points xi = iζ.
La question devient plus réaliste : connaissant plusieurs couples (uk (iζ), f k (iζ)) de
flèches et de forces aux nœuds, peut-on en déduire les valeurs de p(x) et de q(x) ?
Dans la suite, nous noterons cki la valeur approchée de uk (x) au nœud xi = iζ et bki
l’intégrale (II.19) calculée pour f = f k .
Ack = bk . (II.21)
L’espace des matrices carrées Mn×n est euclidien pour le produit scalaire
AH = B (II.24)
Remarque II.3.1
Pour faciliter les calculs, nous chercherons à déterminer les approximations des fonc-
tions p et q aux points xi = iζ, i = 1, ..., n où u et f ont été mesurées et nous
conserverons les fonctions de bases φi définies précédemment.
46 II. Application à l’Identification des Coefficients d’une Equation Différentielle Linéaire
Nous allons donc chercher une approximation des fonctions p et q dans l’espace Vζ :
n
X n
X
pζ (x) = Pk φk (x), qζ (x) = Qk φk (x) (II.25)
k=1 k=1
R
l [p (x)φ0 (x)φ0 + q (x)φ (x)φ (x)]dx = A , 1 ≤ i ≤ n
0 ζ i i ζ i i ii
(II.26)
R l [pζ (x)φ0 (x)φ0 + qζ (x)φi (x)φi+1 (x)]dx = Aii+1 , 1 ≤ i ≤ n − 1
0 i i+1
1 ζ
2ζ
(Pi−1 + 2Pi + Pi+1 ) + 12
(Qi−1 + 6Qi + Qi+1 ) = Aii
(II.27)
− 1 (Pi + Pi+1 ) + ζ
2ζ 12
(Qi + Qi+1 ) = Aii+1
3 1 0 ··· ··· ··· ··· ···
1 6 1 0 ··· ··· ··· ···
0 1 6 1 ··· ··· ··· ···
.. .. .. .. .. .. .. ..
. . . . . . . .
N = .. .. .. .. .. .. .. ..
. . . . . . .
.
··· ··· ··· 0 1 6 1 0
··· ··· ··· ··· 0 1 6 1
··· ··· ··· ··· ··· 0 1 3
1 1 0 ··· ··· ··· ··· ···
0 1 1 0 ··· ··· ··· ···
0 0 1 1 ··· ··· ··· ···
.. ... . . . . . . . . . . . . . . . ..
. .
R=
.. ... . . . . . . . . . . . . . . . ..
. .
··· ··· ··· 0 0 1 1 0
··· ··· ··· ··· 0 0 1 1
··· ··· ··· ··· ··· 0 0 0
et les vecteurs
v = [P1 , ..., Pn ]T , w = [Q1 , ..., Qn ]T
y = [A11 , A22 , ..., Ann ]T , z = [A12 , A23 , ..., An−1n , 0]T ,
le système (II.27) peut devient
M v + ζ 2 N w = 2ζy
6
2
(II.28)
ζ
Rv − Rw = −2ζz.
6
Une condition nécessaire pour que ce système admette des solutions v et w est
l’existence d’un vecteur s = [s1 , ..., sn ]T de Rn tel que
Rs = z.
Il se trouve qu’on peut toujours déterminer une solution particulière s = [s1 , ..., sn ]T
de Rs = z. Le noyau de R étant
−1
1
−1
KerR = R . ,
..
−1
1
M r0 = 0n .
Remarque II.3.2
La précision d’informations à avoir sur les fonctions p et q est liée au nombre de
morceaux dont le matériau est découpé. La matrice A déterminée contient 2n + 3
coefficients ; On connaît la valeur de q à l’origine, donc il faut trouver 2n + 3 valeurs
pour identifier les fonctions p et q.
Tab. II.1 – Les valeurs réelles de la fonction p en différents points et les résultats
obtenus par l’algorithme de projection
II.5 Conclusion
Sommaire
III.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 51
III.2 Position d’un solide dans l’espace . . . . . . . . . . . . . . 53
III.2.1 Vecteur vitesse et vecteur accélération . . . . . . . . . . . 55
III.2.2 Principes fondamentaux de la dynamique . . . . . . . . . 55
III.2.3 Calcul de la résultante dynamique χ . . . . . . . . . . . . 57
III.2.4 Calcul du moment dynamique δ(A) . . . . . . . . . . . . . 57
III.3 δ(A) et χ dans une matrice antisymétrique 4 × 4 . . . . . 58
III.4 Formulation matricielle des équations . . . . . . . . . . . 61
III.4.1 Travail virtuel des efforts . . . . . . . . . . . . . . . . . . 61
III.4.2 Travail virtuel des accélérations . . . . . . . . . . . . . . . 62
III.4.3 Travail virtuel des liaisons "rotation" . . . . . . . . . . . . 63
III.4.4 Equation matricielle du mouvement du solide rigide . . . 63
III.5 Identification de la matrice J . . . . . . . . . . . . . . . . 66
III.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
III.1 Introduction
corps humain.
L’identification des caractéristiques d’inertie d’un athlète peut contribuer à l’amé-
lioration de ses performances, or ces caractéristiques sont généralement mal connues.
Auparavant, les caractéristiques d’inertie des segments du corps humain étaient ob-
tenues par régression linéaire à partir des paramètres inertiels d’une certaine popu-
lation.
On décompose l’athlète en un système composé de 16 troncs de cônes rigides
reliés entre eux par des articulations. L’identification des caratéristiques d’inertie de
l’athlète nécessite donc l’identification de 160 paramètres, un seul solide rigide étant
caractérisé par 10 paramètres.
Afin d’identifier les caractéristiques d’inertie d’un système composé de plusieurs
solides rigides, nous commençons par identifier celles d’un seul solide rigide. Ensuite
nous identifierons celles d’un système de deux solides rigides en tenant compte des
liaisons.
De nombreux auteurs ont identifié les dix caractéristiques d’inertie d’un solide
rigide en les regroupant dans un même vecteur colonne 10 × 1. Ils utilisent ensuite
pour déterminer ce vecteur 10 × 1, le formalisme récursif de Newton-Euler.
Nous avons cherché à améliorer les techniques existantes en commençant par for-
muler matriciellement les équations du mouvement. Dans cette formulation, les 10
caractéristiques d’inertie du solide sont regroupées dans une matrice symétrique 4×4
où la matrice du tenseur d’inertie de Poinsot occupe un bloc 3 × 3. Cette matrice fut
définie antérieurement par Fayet [42] qui la nomma " la matrice du tenseur d’iner-
tie global" ; pour notre part, nous la nommerons "matrice de Poinsot généralisée".
Une première formulation matricielle des équations du mouvement fut introduite en
1996 par Legnani et al.[72]. Dans leur formulation, la matrice de Poinsot généralisée
n’est pas constante, elle varie suivant la position du solide. Nous faisons le choix
d’exprimer la matrice de Poinsot dans le repère lié au solide, ainsi ses composantes
sont constantes et nous pouvons appliquer plus facilement notre méthode d’identifi-
cation. Le mouvement est dans R4 en ajoutant une quatrième coordonnée fixe aux
trois coordonnées d’un point du solide.
Les équations du mouvement sont alors formulées dans l’espace des matrices 4 × 4
antisymétriques.
Pour identifier la matrice 4×4 généralisée de Poinsot, nous avons utilisé l’algorithme
de projection exposé au chapitre I. La méthode a été validée à partir de plusieurs
mouvements simulés. Les simulations du mouvement du solide rigide ont été bruitées
III.2. POSITION D’UN SOLIDE DANS L’ESPACE 53
Considérons un système mécanique composé d’un seul solide rigide. Il est repéré
dans un référentiel R0 (O ; x0 , y0 , z0 ). Le vecteur unitaire k0 est choisi vertical
ascendant et la base B0 (x0 , y0 , z0 ) est choisie orthonormée directe. Nous désignons
par S0 le solide en sa position initiale et par S le solide en sa position à l’instant
t. Soit M0 un point générique du solide S0 . À tout instant t, ce point occupe une
position M dans S. Soit A un point particulier de S qui se trouve en A0 à l’instant
initial. Un second référentiel d’origine A, R(A ; x, y, z), de base orthonormée directe
B(x, y, z), est lié au solide S lors de son mouvement. Les référentiels considérés sont
schématisé sur la figure III.1. Les seules transformations pouvant agir sur le solide
sans le déformer sont les rotations et les translations. Pour passer de sa position
initiale S0 à la position courante S, le solide subit une rotation composée avec une
translation. A l’instant t, nous noterons
T1
→
− −→
T (t) = T2 = OA (III.1)
T3
le vecteur translation et
R11 R12 R13
R(t) = R21 R21 R23 (III.2)
La relation qui définie à tout instant la position d’un point quelconque M du solide
à partir de sa position initiale M0 est donc :
−−→ →
− −−−→
OM (t) = T (t) + R(t)A0 M0 . (III.4)
54 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
en omettant la variable t.
Une matrice rotation R est une matrice qui possède les deux propriétés fonda-
mentales suivantes :
(i) RT R = RRT = I3
(ii) détR = 1.
I3 , RT et dét désignant respectivement la matrice unité de M3 , la matrice transposée
de R et le déterminant.
De la propriété (i) découle la proposition suivante.
Proposition III.2.1
L’inverse d’une matrice étant unique, la matrice transposée d’une matrice rotation
est son inverse :
RT = R−1 .
Les trois colonnes de la matrice rotation R sont les composantes des vecteurs
unitaires x, y, z de la base B exprimées dans la base B0 [35] :
hx, x0 i hy, x0 i hz, x0 i
R = hx, y0 i hy, y0 i hz, y0 i (III.6)
Les produits scalaires apparaissant dans les composantes de la matrice rotation s’in-
terpètent comme les cosinus directeurs des vecteurs de la base B par rapport à la
base B0 . La matrice rotation, appelée dans ces circonstances matrice des cosinus
−−→
directeurs, permet de calculer les composantes du vecteur AM à partir des compo-
−−−→
santes du vecteur A0 M0 .
La matrice rotation R fait passer de la base orthonormée directe du référentiel fixe
B0 (x0 , y0 , z0 ) à la base orthonormée directe B(x, y, z) du référentiel attaché au solide
S.
Soit u un vecteur de R3 ; nous noterons j(u), l’application linéaire antisymétrique,
qui à tout vecteur v de R3 , associe le produit vectoriel
j(u)v = u ∧ v, u ∈ R3 .
III.2. POSITION D’UN SOLIDE DANS L’ESPACE 55
u1 v2 − u2 v1 −u2 u1 0
Le moment dynamique peut être aussi regardé comme une matrice 3 × 3 antisymé-
triques. Ceci est légitime puisqu’une matrice 3 × 3 antisymétrique a 3 composantes
indépendantes comme un vecteur.
En effet si y, v, et w sont trois vecteurs, nous pouvons remarquer que la formule du
double-produit vectoriel
implique
v ⊗ w − w ⊗ v = j(w ∧ v), (III.10)
56 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
= j (δ(A))
Notons F(M ) la densité massique des forces extérieures agissant sur le solide
S autres que la pesanteur et −gz0 la densité massique de la pesanteur (g désigne
l’intensité de la pesanteur) ; alors la résultante des forces extérieures et le moment
dynamique Mex (A) résultant en A sont tels que :
Z
Fex =
(F(M ) − gz0 ) dm
S
(III.12)
−−→ −−→
Z h i
j (Mex (A)) = (F(M ) − gz0 ) ⊗ AM − AM ⊗ (F(M ) − gz0 ) dm
S
ainsi
−−−→
χ = mT̈ + R̈ mA0 G0 (III.16)
En utilisant
−−−→de nouveau l’équation (III.15) et la propriété du produit tensoriel :
−−−→ T
T̈ ⊗ RA0 M0 = T̈ ⊗ A0 M0 R (ou RT est la transposée de R), nous pouvons
reécrire l’équation (III.17) sous la forme :
Z
−−−→ T −−−→ −−−→ −−−→
j (δ(A)) = T̈ ⊗ mA0 G0 R − R mA0 G0 ⊗ T̈ + R̈ A0 M0 ⊗ A0 M0 dm RT
S0
Z
−−−→ −−−→
−R A0 M0 ⊗ A0 M0 dm R̈T ,
S0
(III.18)
R −−−→ −−−→
Nous reconnaissons dans l’intégrale S0 A0 M0 ⊗ A0 M0 dm le tenseur d’inertie de
Poinsot en A0 du solide S dans sa position initiale ; nous le notons J 0 , ainsi :
−−−→ −−−→
j (δ(A)) = T̈ ⊗ mA0 G0 RT − R mA0 G0 ⊗ T̈ + R̈J 0 RT − RJ 0 R̈T . (III.19)
.
58 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
Remarque III.2.1
La résultante dynamique χ calculée par la formule(III.16) est une fonction linéaire
−−−→
des caractéristiques d’inertie m et mA0 G0 du solide. De même, le moment dyna-
mique δ(A) calculé par la formule (III.19) est une fonction linéaire des caractéris-
−−−→
tiques d’inertie mA0 G0 et J 0 .
Remarque III.2.2
(i) Le tenseur d’inertie de Poinsot défini par
−−−→ −−−→
Z
0
J = A0 M0 ⊗ A0 M0 dm
S0
−→
Z
I= −[j(0M )]2 dm.
S
trI
J0 = I − I. (III.20)
2
Si M(A) est le moment des forces extérieures autre que la pesanteur, alors
(
Fex =F − mgz0
−→ −→ (III.22)
j (Mex (A)) =j (M(A)) − gz0 ⊗ mAG + mAG ⊗ gz0
III.3. δ(A) ET χ DANS UNE MATRICE ANTISYMÉTRIQUE 4 × 4 59
sont obtenues des efforts autres que la pesanteur agissant sur le solide au point A.
Des matrices 4 × 4 Λ et Θ contenant les données cinématiques du solide et de
la matrice antisymétrique Υ(A) regroupant la résultante et le moment résultant des
efforts extérieurs, nous formulons matriciellement l’équation régissant à tout instant
t le mouvement du solide :
ΘJ ΛT − ΛJ ΘT = Υ. (III.25)
60 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
Elle correspond au tenseur d’inertie global défini par Fayet [42] pour un solide rigide
isolé.
Cette formulation des équations du mouvement écrite sous forme de produit et
somme de matrices 4 × 4 est linéaire en J . On espère donc que la matrice J sera
facilement identifiable.
Remarque III.3.1
(i) Les six équations du mouvement sont exprimées en termes d’égalité entre matrice
4 × 4 antisymétrique. Une matrice 4 × 4 antisymétrique a six composantes indépen-
dantes, ceci montre que nous avons traduit exactement les équations du mouveme-
ment.
(i) La matrice Υ rassemble des forces et des moments, ses coefficients ne sont
pas homogènes entre eux. De même la matrice J rassemble des moments d’inertie,
des masses et des vecteursproduits d’une masse par une longueur, ses coefficients ne
sont pas homoène non plus. Soit l une longueur de l’ordre de la dimension du solide.
Les matrices Υ, J et Θ peuvent être rendues homogènes respectivement à des mo-
ments, à des moments d’inertie et des inverses de temps au carré :
" #
j (M(A)) lF
Υ(A) =
−lFT 0
" −−−→#
J0 mlA0 G0
J = −−−→
mlA0 G0 T ml2
" #
R̈ T̈ + gz0 /l
Θ= .
0T 0
Soient F(M ) la densité des forces extérieures autres que la pesanteur et g l’in-
tensité de la pesanteur. Le travail virtuel des efforts extérieurs qui agissent sur le
R
solide S : τ = S (F (M ) − gz0 ) .δM peut se transcrire :
Z " # "−−→#!
F(M ) − gz0 OM
τ= .δ dm
S 0 l
62 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
Le Principe des travaux virtuels impose que : " Le travail virtuel des quantités
d’accélérations est égal au travail virtuel des efforts extérieurs pour tout déplacement
virtuel rigidifiant. Le caractère rigidifiant du déplacement virtuel peut se résumer
par τ ” = 0. D’après la technique des multiplicateurs de Lagrange,
τ0 = τ + τ”
Retranchons gz0 aux forces extérieurs revient au même que de rajouter gz0 aux
accélérations selon l’expérience de pensée d’Einstein dite de l’ascenseur en chute
libre. " #
gz0 /l
Si on ajoute à l’accélération généralisée (III.28), on obtient
0
" # " #
a(M ) + gz0 /l R̈ T̈/l + gz0 /l
= T .
0 0 0
" # "−−→# "−−→# " #!
R F(M ) AM AM F(M )
La matrice S ⊗ − ⊗ dm
0 l l 0
se développe successivement en
"R h −−→ −−→ i R #
S
F(M ) ⊗ AM − AM ⊗ F(M ) dm l S
F(M )dm
R T
−l S F(M )dm 0
" R −−→ R #
j S AM ∧ F(M )dm l S F(M )dm
= R T
−l S F(M )dm 0
" R #
j (M(A)) l S F(M )dm
= R T . (III.30)
−l S F(M )dm 0
III.4. FORMULATION MATRICIELLE DES ÉQUATIONS 65
Proposition III.4.1
Le tenseur de Poinsot généralisé J est une matrice défini-positive
−−−→
Z 2
hV, J Vi = hA0 M0 , Ui + lv dm
S0
Remarque III.4.1
La positivité de la matrice J généralise la positivité du tenseur d’inertie de Poinsot
J 0.
66 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
Considérons le problème :
Á partir des matrices 4 × 4 Λ et Θ obtenues par l’enregistrement du mouvement
d’un solide à différents instants, et de la matrice antisymétrique Υ composée des
forces et moments externes, trouvons la matrice 4 × 4 symétrique et défini-positive
J telle que :
ΘJ ΛT − ΛJ ΘT = Υ. (III.31)
10 X
X 10 10
T T
X
hi tr Sj ΘSi Λ − ΛSi Θ = tr [Sj Υ] . (III.32)
j=1 i=1 j=1
La matrice de composantes Sj ΘSi ΛT − ΛSi ΘT est défini-positive.
Différents types de mouvements ont été simulés à l’aide du logiciel Solid Dyna-
mics. Les simulations représentaient des cas simples : un parrallélépipède de masse
m = 0, 8kg et de côtés 0, 4m, 0, 3m et 0, 5m est soumis à une force ponctuelle agis-
sant en f en A (voir figure III.2). Les composantes des forces et des déplacements
enregistrés sont données par la figure III.3. Ainsi,
" à chaque
# instant et pour chaque
[0] Fex
simulation, la matrice Υ(A) est de la forme : . Seule la valeur de de la
−FT ex 0
force ponctuelle f en A varie d’une simulation à l’autre.
masse 0, 8 kg
0, 20
−→
position du centre de masse AG = 0, 15 m
0, 25
0, 02267 0 0
matrice d’inertie en G I(G) = 0 0, 02733 0 m2 kg
0 0 0, 01667
Cinq mouvements sont simulés : pour deux mouvements, les forces extérieures
agissant au point A sont restées relativement faible, et pour trois autres mouvements
simulés, les variations des forces extérieures sont importantes.
Les résultats sont dans les tableaux suivants. Les analyses portent sur les résultats
obtenus des données bruitées.
Nous avons examiné les rapports kJ − Jid k/kJ k afin d’apprécier les résultats,
Jid étant la matrice du tenseur d’inertie identifiée. Pour :
68 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
Tab. III.1 – Résultats de la simulation avec une force maximale de 91,6N ; la norme
de l’écart entre la matrice identifiée et la matrice calculée J est 0, 0708
Tab. III.2 – Résultats de la simulation avec une force maximale de 49.26N ; la norme
de l’écart entre la matrice identifiée et la matrice calculée J est 0.0782
Tab. III.3 – Résultats de la simulation avec une force maximale de 21.15N ; la norme
de l’écart entre la matrice identifiée et la matrice calculée J est 0.2893
Tab. III.4 – Résultats de la simulation avec une force maximale de 57.3N ; la norme
de l’écart entre la matrice identifiée et la matrice calculée J est 0.1022
III.5. IDENTIFICATION DE LA MATRICE J 69
Tab. III.5 – Résultats de la simulation avec une force maximale de 23.6N ; la norme
de l’écart entre la matrice identifiée et la matrice calculée J est 0.2018
Tab. III.6 – Résultats obtenus par explotation simultanée des données bruitées des 5
simulations et par explotation simultanée des données bruitées des forces de 49,26N
et de 91,6N
70 III. Identification des Caractéristiques d’Inertie d’un Solide Rigide
III.6 Conclusion
Deuxième partie
Sommaire
IV.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 73
IV.2 Problème direct . . . . . . . . . . . . . . . . . . . . . . . . 75
IV.2.1 Généralités et définitions . . . . . . . . . . . . . . . . . . . 75
IV.2.2 Equation de la chaleur dans un matériau . . . . . . . . . . 78
IV.2.3 Condition initiale et conditions aux limites . . . . . . . . . 80
IV.2.4 Résolution du problème direct par la méthode des éléments
finis (MEF) . . . . . . . . . . . . . . . . . . . . . . . . . . 84
IV.3 Problème inverse . . . . . . . . . . . . . . . . . . . . . . . . 88
IV.3.1 Techniques de mesures de champs de température . . . . . 90
IV.3.2 Méthodes de résolution des problèmes inverses . . . . . . 91
IV.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
IV.1 Introduction
Avant d’aborder le problème inverse, nous commençons par une étude prélimi-
naire du problème direct. Au cours de cette étude, nous décrirons le type de transfert
considéré dans nos travaux et définirons les paramètres thermophysiques caracté-
risant les matériaux homogènes et isotropes. Nous introduirons quelques notions
essentielles relatives à la conduction de la chaleur dans un solide opaque.
74 IV. Introduction aux Problèmes Direct et Inverse en Thermique
∂2T 1 ∂T
2
=
∂x a ∂t
Avant d’aborder les généralités, donnons quelques définitions des notions fonda-
mentales [52] :
(i) Variable d’état : L’ensemble des propriétés physiques que l’on peut repré-
senter par une série de grandeurs macroscopiques α (volume, masse, température,
déformations...) caractérisent l’état d’un système.
(ii) Fonction d’état : Une fonction d’état décrit les propriétés du système et
traduit par des égalités ou inégalités des lois de la thermodynamique.
(iii) Température : La température T est une variable d’état permettant de
définir la notion d’équilibre thermique. C’est une fonction scalaire des coordonnées
des points du domaine et du temps t.
(iv) Capacité thermique : La capacité thermique cα d’un corps est une gran-
deur permettant de quantifier la possibilité qu’a un corps d’absorber ou de restituer
de l’énergie par échange thermique au cours d’une transformation pendant laquelle
sa température varie.
(v) Flux thermique et densité de flux thermique : Considérons un matériau
parallélépipédique occupant un domaine Ω et chauffé de manière uniforme sur une
de ses surfaces Σ. Soit Q la quantité de chaleur reçue par un élément de surface dΣ.
On définit [86] :
Cette loi qui traduit la relation linéaire entre la densité de flux thermique et le
gradient de température est dite Loi de Fourier (1882).
L’opérateur k̄¯, qui définit la conductivité thermique du milieu considéré, est un
tenseur du second ordre dépendant de la température T du point M . Par exemple,
pour une conduction anisotrope comme dans les matériaux fibreux, k n’a pas la
même valeur dans toutes les directions. Pour un milieu isotrope et homogène, la
conductivité thermique se caractérise simplement par un scalaire k dépendant de la
→
− −−→
température T : ψ = −k grad(T ). De plus, si les écarts de température ne sont pas
trop élevés alors on peut considérer la conductivité k constante pour le milieu.
→
−
Par convention [86], ψ est compté positivement dans le sens d’écoulement de la
chaleur, donc dans le sens des températures décroissantes.
Pour un milieu unidimensionnel x, la loi de Fourier est traduite par :
dT
ψ = −k où k est la conductivité thermique.
dx
(vii) Energie interne : L’énergie interne E d’un système de domaine d’extension
Ω est reliée à son énergie interne spécifique e par la relation :
Z
E= ρedΩ.
Ω
IV.2. PROBLÈME DIRECT 77
Z
S= ρsdΩ (IV.4)
Ω
−
→
Z Z
dΨ ds dT
div ψ dΩ = ρ +T + s dΩ
Ω Ω dt dt dt
∂ 2 Ψ dT
Z
dΨ dT ∂Ψ
= ρ −T − dΩ
Ω dt ∂T 2 dt dt ∂T
∂ 2 Ψ dT
Z
∂Ψ dαj dT ∂Ψ
= ρ −T − dΩ, j = 1, ..., n
Ω ∂αj dt ∂T 2 dt dt ∂T
αj , j = 1, ..., n sont les variables d’état définis plus haut.
La capacité calorifique c définie à αj constant est
∂2Ψ
c = −T .
∂T 2
dT /dt représente la dérivée particulaire de la température. C’est la dérivée suivant
le mouvement de déformation et de l’évolution temporelle du matériau. Elle s’écrit :
dT −−→ ∂T
=−
→
u .grad(T ) + ,
dt ∂t
où −
→u est la vitesse des volumes élémentaires du matériau en déformation.
Certaines études [89] ont montré que dans le cas de solicitations uniaxiales homo-
gènes à faible vitesse de déformation dans les matériaux malgré la faible conductivité
−−→
thermique, par exemple comme les polymères, le terme − →u .grad(T ) est négligeable
devant ∂T /∂t.
De même pour des variations très faible de température, les déformations et les effets
sur la microstructure (énergie bloquée) sont négligeables.
De ces hypothèses réductrices, le bilan d’équilibre devient
→
−
Z Z
∂T
div ψ dΩ = ρc dΩ.
Ω Ω ∂t
Cette équation aux dérivées partielles linéaires, du deuxième ordre, admet mathéma-
tiquement une infinité de solutions faisant intervenir des constantes arbitraires. Elle
n’a de sens physique que pour des conditions, appliquées à un domaine d’espace-
temps. Ces conditions sont ici les causes qui déterminent l’évolution du phénomène.
Pour déterminer l’évolution physique de la température dans Ω, il est indispen-
sable de connaître la répartition initiale de la température en tout point du domaine
et aussi la loi de variation de la température T et/ou de sa dérivée normale sur
la frontière du domaine. Ces conditions interprétées par des équations sont dites
condition initiale et conditions aux limites.
Ce sont des conditions que l’on impose sur la surface frontière ou limite Σ. Elles
traduisent mathématiquement les conditions physiques réelles imposées. En ther-
mique, il existe différents types de conditions aux limites. Les plus fréquentes ren-
contrées sont les suivantes :
Température imposée
Ce sont des problèmes dits de Dirichlet, où la température est une fonction donnée
d’un point M et du temps t. Elle est dite condition d’isothermie et réalisée par
contacte du domaine avec des thermostats :
TΣ = f (MΣ , t) en tout point MΣ de la frontière de Σ.
Flux de chaleur imposé
La densité de flux thermique est une fonction donnée du point à la frontière et/ou
du temps : ces problèmes sont dits de Neumann et
− −
→ −−→
ψ .→
n = −k grad(T ).−
→
n = f (MΣ , t) en tout point MΣ de la frontière Σ. (IV.11)
Par analogie avec la loi de Newton, hr est appelé coefficient d’échange thermique par
rayonnement.
En utilisant la loi de Fourier et la continuité de flux à l’interface d’échange, la
condition de transfert par rayonnent s’écrit :
−−→
−k grad(T ).−
→
n = σF (TΣ4 − Ta4 ) ≈ hr (TΣ − Ta ) (IV.15)
Σ1 ∪ Σ2 = Σ et Σ1 ∩ Σ2 = ∅
−−→
Imposons alors sur la frontière Σ1 la condition de flux imposé : −k grad(T ).−→n =q
−−→ →
−
et sur Σ2 une condition de tranfert par convection : −k grad(T ). n = h(T − Ta ).
Dans ce cas, si on adjoint la condition initiale à ces conditions, on peut décrire un
phénomène de transfert de chaleur dans un matériau par les équations
−−→
pour (− →
ρc ∂T
∂t
− k div(grad(T )) = 0 x , t) ∈ Ω × R+
−k − −→
grad(T ).− →
n =q en tout point de Σ1 , t > 0
−− → (IV.16)
−k grad(T ). →
−
n = h(T − T ) en tout point de Σ t > 0
a 2
−
T (→x , 0) = T0 (−→ pour − →
x) x ∈Ω
Remarque IV.2.2
Nous avons choisi ces conditions aux limites plutôt que celles de Dirichlet (T (x, t)
connues sur Σ2 ) car nous voulons faire apparaître les termes de convection h.
Le problème (IV.16) est donc constitué d’une équation aux dérivées partielles
munies de conditions aux limites et d’une condition initiale. A cause de la présence
de conditions aux limites, on dit que le problème (IV.16) est un problème aux limites ;
mais on dit aussi que c’est un problème de Cauchy à cause de la donnée initiale en
temps.
Dans le cadre de la résolution du problème (IV.16), on peut y procéder directement
pour déterminer le champ de trempérature T (− →
x , t). Ceci peut se faire de façon
analytique, ou en utilisant la transformation de Laplace. Néanmoins, la plupart du
temps on ne peut en déterminer qu’une solution approchée. Il existe de nos jours
des méthodes numériques donnant plus d’informations sur le champ de température
84 IV. Introduction aux Problèmes Direct et Inverse en Thermique
T (−
→
x , t). Parmi ces méthodes numériques, la méthode des éléments finis est la plus
utilisée.
La MEF est une des méthodes les plus efficaces et les plus populaires pour ré-
soudre numériquement des problèmes aux limites. Le principe de cette méthode
peut être basé sur l’approximation variationnelle. La formulation variationnelle
est posée dans un espace de Hilbert (dans ce document c’est l’espace de Sobolev
H 1 (Ω) [58]). Rappelons que les espaces de Sobolev qui sont des espaces "naturels"
de fonctions permettant de résoudre les formulations variationnelles d’équations aux
dérivées partielles. Le lecteur intéressé, ou tout simplement curieux, trouvera les
idées et les arguments clés qui lui permettront de bien comprendre la structure et
l’intérêt des espaces de Sobolev dans le second chapitre ou plus en détail dans les
ouvrages [58] et [25].
L’idée de base est de remplacer l’espace H 1 (Ω) de dimension infinie par un sous-
espace Vζ de dimension finie de H 1 (Ω). Il faut d’une part construire une base φi des
fonctions à support borné, c’est-à-dire des fonctions non nulles uniquement sur une
petite partie du domaine de définition, et d’autre part, la réunion des supports des
fonctions φi doit recouvrir Ω. Ces fonctions sont généralement appelées fonctions
"test".
Remarque IV.2.3
Dans la pratique, il faut que la réunion des supports des fonctions φi soit la plus
proche possible du domaine Ω. Cependant, il est difficile voire impossible dans le
cas des figures à géométrie très complexe de reconstituer tout le domaine Ω par la
réunion des supports des fonctions test.
Cette difficulté s’élimine dans les cas unidimensionnels et la méthode des éléments
finis conduit à des solutions exactes [30].
La formulation variationnelle consiste à trouver une forme bilinéaire a(., .), une forme
linéaire L(.), et un espace de Hilbert V = Vζ tels que le problème de conduction de
la chaleur soit équivalent à :
T : ]0, T [ −→ H 1 (Ω)
t −→ T (t),
La formulation de Galerkin est définie par deux intégrales représentant les poids
résiduels [3] et [70]. La méthode résiduelle pondérée se traduit en choisissant un
élément φj de Vζ par :
−−→ −−→
Z Z
∂T
φj k grad(T ).−
→
φj k div(grad(T )) − ρc dΩ − n + q dΣ = 0 (IV.18)
Ω ∂t Σ
Ces équations sont généralement écrites sur chaque domaine élémentaire, et après
les intégrations, on passe à l’assemblage de chaque élément pour obtenir un système
linéaire [58] et [67]. Nous définirons des fonctions élémentaires φi (x) permettant
d’avoir directement le système linéaire attendu.
Parmis les équations de conduction que nous utiliserons, une condition aux li-
mites de Dirichlet sera imposée. Les conditions aux limites de type Dirichlet peuvent
correspondre à la situation où des phénomènes physiques différents agissent sur des
portions distinctes de la frontière. Si nous reprenons Σ2 et Σ1 comme définis précé-
demment, et nous imposons les conditions
−−→
−k grad(T ).−
→
n =q sur Σ1
−−→
k grad(T ).−
→
n = h(T − Ta ) sur Σ2 ,
C θ̇ + Kθ = F (IV.21)
où
– (i) θ̇ désigne le vecteur dérivé temporaire du vecteur température θ,
– (ii) C désigne la matrice de capacité, aussi appelée matrice d’amortissement,
– (iii) K est la matrice de conductance, encore appelée matrice de raideur,
– (iv) F est le vecteur flux résultant des conditions aux limites.
Les composantes des matrices C, K et du vecteur F sont données par
Z
Cij = ρcφj φi dΩ (IV.22)
Ω
−−→ −−→ −−→
Z Z Z
Kij = →
−
k grad(φj ).grad(φi )dΩ − φj k grad(φi ). n dΣ + φj hφi dΣ (IV.23)
Ω Σ Σ2
Z Z
Fj = − φj qdΣ + φj hTa dΣ (IV.24)
Σ1 Σ2
Remarque IV.2.4
La matrice C est symétrique et défini-positive ; la matrice K est symétrique et posi-
tive. Les matrices C et K sont tridiagonales.
θ(0) = [θ1 (0) θ2 (0) ... θn (0)] = [T0 (0) T0 (ζ) ... T0 (iζ) ... T0 (L)]
Les travaux de Bell ont été repris en 1961 par Parker et al. [103] en mettant
au point la méthode flash. La méthode flash dans son principe, est une technique
impulsionnelle qui consiste à soumettre la face avant d’un échantillon plan à une
impulsion de flux de chaleur de courte durée et à observer l’évolution temporelle de
la température en un ou plusieurs points de la face arrière de l’échantillon. Outre
l’identification de la diffusivité thermique, elle permet d’accéder indirectement à la
conductivité thermique lorsque la capacité massique et la masse volumique sont
connues.
Les méthodes d’estimation de paramètres liées à la méthode flash ont fait l’objet
de nombreuses études [11] et [100]. Ces méthodes ont été revisitées en tenant compte
des outils liés aux méthodes d’estimation de paramètres (méthodes inverses). Les
trois principaux types de méthodes flash sont présentées ci-après :
92 IV. Introduction aux Problèmes Direct et Inverse en Thermique
La méthode du plan chaud est une extension de la méthode du fil chaud à une
géométrie plane. Elle permet d’identifier l’effusivité thermique du matériau [38].
Le principe est le suivant : un flux de chaleur uniforme est imposé à l’interface de
deux échantillons symétriques de section quelconque et d’extension infinie dans la
direction perpendiculaire à l’interface. Ceci est réalisé par la mise en place d’un élé-
ment chauffant mince occupant toute la section de l’interface. Les échantillons sont
assimilables à un milieu infini si le rapport de leur longueur à l’épaisseur de l’élément
chauffant est supérieur à 20. Par ailleurs, les faces latérales des deux échantillons
sont isolées et le transfert peut être considéré comme unidirectionnel. La réponse
en température au cours du temps est mesurée dans le plan du chauffage par un
thermocouple. Le comportement asymptotique aux temps longs, du thermogramme
ainsi obtenu, est proportionnel à la racine carrée du temps, le coefficient de propor-
tionnalité étant directement lié à l’effusivité des échantillons. Le comportement aux
temps courts est influencé par l’inertie thermique de la sonde (élément chauffant +
thermocouple) et par la résistance de contact sonde-milieu.
Cette méthode permet de mesurer l’effusivité thermique avec une précision de l’ordre
de ±6% pour peu que l’intervalle de temps sur lequel est réalisée l’identification soit
correctement choisi de sorte que puissent être négligés les effets inertiels de la sonde
aux temps courts et les éventuelles pertes thermiques aux temps longs.
94 IV. Introduction aux Problèmes Direct et Inverse en Thermique
IV.4 Conclusion
La méthode des éléments finis qui est basée sur la technique des projections,
permet donc d’adapter le mieux possible le maillage au domaine géométrique. Elle
conduit à une réduction considérable des erreurs de résolution numérique. La for-
mulation de Galerkin nous a emmené à écrire l’équation aux dérivées partielles de
conduction en une équation différentielle du premier ordre en dimension finie. Ces
circonstances expliquent notre engouement manifesté pour l’application de la mé-
thode des éléments finis pour la résolution de l’équation de conduction de la chaleur.
Quelle que soit la méthode utilisée, les différentes techniques relient le modèle
théorique à l’observable physique c’est-à-dire aux paramètres de contrôle comme la
température. Les mesures de température sur un échantillon sont ensuite injectées
dans les équations régissant la conduction afin de remonter aux grandeurs inconnues
recherchées. Toutes ces méthodes conduisent obligatoirement à la minimisation d’une
fonction coût.
95
Chapitre V
Sommaire
V.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 96
V.2 Problème direct . . . . . . . . . . . . . . . . . . . . . . . . 97
V.2.1 Equation de conduction de la chaleur . . . . . . . . . . . . 97
V.2.2 Formulation du problème direct . . . . . . . . . . . . . . . 98
V.3 Résolution du problème direct . . . . . . . . . . . . . . . . 98
V.3.1 Discrétisation en espace par MEF . . . . . . . . . . . . . . 98
V.3.2 Intégration dans le temps . . . . . . . . . . . . . . . . . . 101
V.3.3 Algorithme de résolution du problème direct . . . . . . . . 102
V.4 Choix des pas temporel et spatial . . . . . . . . . . . . . . 103
V.4.1 Stabilité de l’algorithme du problème direct . . . . . . . . 103
V.4.2 Tests de stabilité . . . . . . . . . . . . . . . . . . . . . . . 105
V.5 Formulation du problème inverse . . . . . . . . . . . . . 106
V.6 Procédure d’identification des paramètres thermophy-
siques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
V.6.1 Méthode des moindres carrés . . . . . . . . . . . . . . . . 106
V.6.2 Calcul de précision sur les paramètres identifiés . . . . . . 110
V.6.3 Construction d’un algorithme de résolution du problème
inverse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
V.7 Validation de l’algorithme . . . . . . . . . . . . . . . . . . 113
V.7.1 Simulation du champ de température . . . . . . . . . . . . 114
V.7.2 Exploitation du champ de température simulé . . . . . . . 115
V.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
V.1 Introduction
∂T ∂2T
ρc =k 2, x ∈ [0, L], t ≥ 0. (V.1)
∂t ∂x
À cette équation, on adjoint les conditions aux limites de flux
∂T
−k = q, x=0 t>0 (V.2a)
∂x
∂T
k = 0, x=L t>0 (V.2b)
∂x
et la condition initiale
Remarque V.3.1
Les fonctions test ainsi définies sont bien dérivables au sens des distributions et leurs
dérivées sont continues. Donc les fonctions φi sont des éléments de Vζ .
Fi = qφi (0).
Remarque V.3.2
La somme des composantes de la matrice K suivant chaque ligne est nulle. On peut
en conclure que le noyau de la matrice K est composé des vecteurs colinéaires au
vecteur de composantes toutes égales à 1 ([1 1 ... 1 1]T ) de Rn+1 .
Remarque V.3.3
On peut obtenir l’équation (V.7) par la formulation variationnelle et dan ce cas.
Les calculs des composantes des matrices C et K et du vecteur F sont détaillés dans
l’annexe A.
Remarque V.3.4
Si les matrices C et K et le vecteur F dépendent de la variable t du temps, il faut les
calculer aussi à chaque instant.
Le vecteur dérivé température θ̇(t) est obtenu directement par l’équation différen-
tielle (V.7).
102 V. Construction d’un Algorithme d’Identification de Paramètres Thermophysiques
et
N = C + ν∆tK
D 0 = R0 .
4. calculer le gradient Ri+1 = −S(θi+1 ), puis évaluer hRi+1 , Ri+1 i. S’il n’est pas
nul, continuer
5. calculer
h(Ri+1 − Ri ), Ri i
γi =
hRi , Ri i
V.4. CHOIX DES PAS TEMPOREL ET SPATIAL 103
Remarque V.3.5
La matrice N étant défini-positive, l’algorithme converge au plus en (n+1)-itérations
[84].
Le vecteur θ̇(t) = C −1 (F − Kθ(t)) est calculé avant le calcul de S(θ). Autrement dit,
chaque fois qu’une solution θ(t) est déterminée, il faut calculer θ̇(t) avec (V.5) afin
de pouvoir l’incorporer dans (V.8) pour la détermination de θ(t + ∆t).
∂lf
= fu...u f est indicée en u l fois
∂ul
De la formule de Taylor avec le reste de Lagrange appliquée à T (x, t) pour (1 ≤ i ≤
n − 1) et (1 ≤ j ≤ n − 1), on a :
? d’une part à l’ordre 2 en t
1
T (xi , tj+1 ) = T (xi , tj ) + Tt (xi , tj )∆t + Ttt (xi , τ )(∆t)2 , τ ∈ [tj , tj+1 ] (V.9)
2
ou encore
T (xi , tj+1 ) − T (xi , tj ) 1
Tt (xi , tj ) = − Ttt (xi , τ )∆t, τ ∈ [tj , tj+1 ]
∆t 2
104 V. Construction d’un Algorithme d’Identification de Paramètres Thermophysiques
R(x, t) = − ρc[Ttt (xi−1 , tj )(∆t) + Ttt (xi , tj )(∆t) + Ttt (xi+1 , tj )(∆t)]+
O(∆t)2 + kTxxxx (xi , tj )ζ 2 + O(ζ 2 )
Il convient donc de choisir le pas de temps ∆t et le pas spatial ζ tels que −6ρc +
kζ 2 < 0 pour réduire le reste R(xi , tj ) de l’équation (V.12) inférieure ou égale à
O(∆t)2 + O(ζ 3 ).
Pour le calcul de dθ/dt dans l’équation (V.7) avec la formule de Crank-Nicolson, il
faut que ∆t et ζ vérifient
6k∆t
1< . (V.13)
ρcζ 2
Si la condition (V.13) n’est pas respectée, alors la solution de l’équation (V.7) pré-
sentera des oscillations spatiales contraires à la physique. Cette condition concerne
généralement l’épaisseur ζ de la première maille concernée par le choc thermique.
Cette contrainte donne la longueur maximale en fonction de l’incrément de temps
utilisé : s
6k∆t
ζ ≤ ζmax = .
ρc
Avec la formule des trapèzes généralisées, l’algorithme d’intégration converge si
s
6νk∆t
ζ ≤ ζmax = . (V.14)
ρc
On retrouve une autre démonstration de ce critère faite par J.-M. Bergheau et al.
[67].
Pour valider le critère de stabilité (V.14), nous avons effectué deux applications
numériques. Nous avons choisi les valeurs fixes suivantes :
– longueur de la plaque : L = 0.04 m,
– pas de discrétisation : ζ = 6.6 10−4 m,
– conductivité thermique : k = 0.45 W/m/◦ C,
– température initiale : T0 = 20◦ C,
– chaleur volumique : ρc = 1.74 106 J/m3 /◦ C
– puissance de chauffage constante : q = 427.355 W/m.
2
Pour ces valeurs choisies, le pas de temps (seuil inférieur) est de ∆ts = ρcζ 6k
≈
0.57 s. Nous avons choisi deux pas de temps d’intégration : dans un premier temps,
un pas de temps supérieur au pas de temps seuil puis un second pas de temps inférieur
au pas de temps seuil. On constate sur les figures (Fig.V.3) que les solutions obtenues
dépendent du pas de temps choisi. Avec ∆t = 1 s, la solution est relativement
106 V. Construction d’un Algorithme d’Identification de Paramètres Thermophysiques
où k .k est la norme associée au produit scalaire h., .i dans Rn+1 , tf est la durée des
mesures.
En utilisant la définition de la norme et la linéarité du produit scalaire, la fonction
PROCÉDURE D’IDENTIFICATION DES PARAMÈTRES
V.6.
THERMOPHYSIQUES 107
La résolution du système (V.17) permet d’estimer les valeurs des paramètres ther-
mophysiques k et ρc.
Remarque V.6.1
On peut avoir directement le système linéaire (V.17) en multipliant l’équation
dθ
ρcC + kKθ(t) = F
dt
à droite par la transposée du vecteur θ(t) puis on intègre sur la durée des mesures et
on obtient
Z tf Z tf Z tf
dθ
ρc hC , θ(t)idt + k hKθ(t), θ(t)idt = hF, θ(t)idt
0 dt 0 0
A11 u1 + A12 u2 = f1 .
A21 u1 + A22 u2 = f2
Au = f .
Remarque V.6.2
L’équation (V.17) est mal conditionnée. En effet,
(i) les composantes A11 , A12 et A22 de la matrice A sont respectivement homo-
gènes à (m.◦ C/s)2 , (◦ C)2 /s et (◦ C/m)2 ,
(ii) u1 et u2 du vecteur u sont respectivement homogènes à J/m3 /◦ C et W/m/◦ C
(iii) f1 et f2 du vecteur f sont respectivement homogènes à W.◦ C/s et W.◦ C/m2 ,
alors qu’une matrice ou un vecteur doit être homogène.
PROCÉDURE D’IDENTIFICATION DES PARAMÈTRES
V.6.
THERMOPHYSIQUES 109
Par changement d’échelle, nous avons rendu homogène la matrice A, les vecteurs u
et f . Le conditionnement donne :
1 dθ dθ 1 dθ
hC i hC
Z tf L2 dt , C dt tf dt
, Kθi Z tf hC dθ
dt
, Fi
A= dt f = dt
0 1 dθ L2 0 L2
tf
hC dt , Kθi t2 hKθ, Kθi tf
hKθ, Fi
f
hv, M vi = 0 ⇒ v = 0.
(ii) On ne peut parler d’erreur que si l’on dispose d’une valeur de référence que l’on
peut considérer comme valeur vraie. Le calcul d’erreur est possible lorsqu’il s’agit de
résultats obtenus par simulation numérique, et dans ce cas, est considérée comme
valeur vraie, toute valeur d’entrée.
Vu que pour la plupart des mesures, exemple des mesures de caractérisation de ma-
tériaux, on ne dispose pas souvent de valeur de référence, par des analyses objectives
des moyens utilisés tout au long des calculs, nous pouvons nous faire une idée de
l’erreur maximale commise : c’est l’incertitude.
Une mesure est complètement indiquée par la valeur estimée et l’intervalle à l’in-
térieur duquel on est à peu près certain que se situe la valeur vraie. L’incertitude
absolue est, par définition, la demi-longueur de cet intervalle. Il existe deux manière
de procéder au calcul d’incertitude absolue.
Lorsque la distribution du bruit est supposée faiblement aléatoire, on calcule des
incertitudes maximales, alors qu’on calcule des incertitudes probables lorsque la dis-
tribution est supposée fortement aléatoire et indépendantes.
Remarque V.6.4
Les incertitudes sont calculées à partir des écarts types sur les données ou mesures.
Les incertitudes probables δprob calculées en prennant 2σ sur les mesures expriment
68% de chance d’être dans l’intervalle de confiance déterminé. Avec une incerti-
tude 2δprob , on a 95% de chance pour que la valeur estimée soit dans l’intervalle de
confiance.
Des hypothèses faites sur la distribution des bruits, il est souhaitable de calculer les
incertitudes probables.
Les composantes de la matrice A et du vecteur f étant obtenues par produit
scalaire de vecteurs, nous calculerons les incertitudes probables sur le réel
Z tf
a= hx, yidt.
0
R tf
Par exemple, la première composante de la matrice A est A11 = 0
hCdθ/dt, Cdθ/dtidt,
et donc x = y = Cdθ/dt.
112 V. Construction d’un Algorithme d’Identification de Paramètres Thermophysiques
m−1
X (tj+1 − tj )
= hxj , yj i + hxj+1 , yj+1 i , avec xj = x(tj )
j=0
2
De la formule (V.18), on calcule les incertitudes probables sur chacunes des compo-
santes de A et f .
Remarque V.6.5
Des hypothèses faites sur les mesures de température, les incertitudes relatives sur
chaque mesure de température en tout point i est δθ i.e. δ(T (xi , t) = δθ. Il vient que
δ(A22 ) = 0 puisque tout vecteur colinéaire au vecteur de composantes toutes égales
à 1 est dans le noyau de la matrice K.
Des égalités
Par exemple, l’incertitude probable δ(k) est ainsi obtenue avec un niveau de confiance
de 95%, c’est-à-dire que nous sommes certains 95 fois sur 100 que notre meilleure
valeur estimée pour la conductivité thermique est comprise dans l’intervalle [k −
δ(k), k + δ(k)]. Il en est de même pour le coefficient de la chaleur volumique.
Finalement les résultats seront présentées de la façon suivante :
Après avoir établi les étapes informatiques de base requises pour la solution du
problème inverse ci-dessus, nous présentons maintenant l’algorithme permettant sa
résolution.
Pour déterminer simultanément les paramètres thermophysiques k et ρc à partir
d’un champ de température, il faut :
1. calculer par intégration numérique les troix coefficients de la matrice A ;
2. calculer par intégration numérique les deux coefficients du vecteur f ;
3. déterminer algébriquement la solution u de l’équation Au = f .
Remarque V.6.7
Les mesures n’étant pas continues, on peut remplacer les intégrales par des sommes.
Nous avons opté pour la méthode des trapèzes pour calculer les intégrales. Il existe
d’autres méthodes d’intégration pour pouvoir calculer les composantes de la matrice
A et du vecteur f .
La matrice A et le vecteur v sont appelés généralement matrice de mesures et
vecteurs de mesures.
ksi = 0, 45 W/m/◦ C
La température initiale est T0 (x) = 20◦ C et le flux de chaleur auquel est soumis la
plaque est choisi constant qs = 397W/m.
Remarque V.7.1
Les moyens expérimentaux dont nous disposons ne nous permettront pas d’évaluer
le flux de chaleur variable dans le temps. Le flux de chaleur sera, dans les condi-
tions expérimentales, la puissance électrique par unité de longueur fournie par un
générateur.
θi = Θi + aω
Les mesures de température ont été effectuées pendant une durée tf . Supposons
enregistrer m mesures de température à intervalle de temps constant ∆t c’est-à-dire
tf = m∆t. Le vecteur θ(t) n’est donc connu qu’aux instants 0, ∆t, ... i∆t, ..., tf .
Nous définissons trois méthodes de dérivation première, et à titre d’exemple pour
a = 0.02◦ C (c’est l’écart-type des mesures effectuées par les caméras infrarouges
actuelles), nous avions analysé les profils des vecteurs dérivés de chaque méthode.
? dérivée première décentrée à droite
Remarque V.7.3
Cette formule ne permet pas de calculer la dérivée de θ en t = 0.
La figure Fig.V.6 illustre l’allure des courbes dθ(t)/dt obtenues à partir de champs
de température bruité et non bruité pour un même instant t. Il ressort des analyses
que la dérivée centrée conviendrait mieux dans les techniques d’intégration et de
dérivation numérique. Pour ne pas perdre des informations, la dérivée à l’instant
t = 0 est recupérée par la dérivée centrée à gauche alors que la dérivée à t = tf est
obtenue par la dérivée centrée à droite. Le graphe "Les trois Dérivées" correspond
aux courbes obtenues par dérivation des données non brutées ; elles sont identiques.
3k∆t
> 1.
ρcζ 2
Fig. V.7 – Ecarts entre valeur simulée et valeurs identifiées pour la chaleur volu-
mique ; nombre d’intervalles et durée de contrôle du flux nécessaires pour une bonne
estimation
V.7. VALIDATION DE L’ALGORITHME 117
Fig. V.8 – Ecarts entre valeurs simulées et valeurs identifiées pour la conductivité
thermique ; nombre d’intervalles et durée de contrôle du flux nécessaires pour une
bonne estimation
Autrement dit, nous sommes certains 95 fois sur 100 que la conductivité thermique et
118 V. Construction d’un Algorithme d’Identification de Paramètres Thermophysiques
Tab. V.1 – Résultats des paramètres estimés à partir d’un champ de température
simulé
la chaleur volumique estimées se situent respectivement dans les intervalles [0, 404 −
0, 472] et [1669000 − 1809000].
V.8 Conclusion
La MEF est utilisée pour écrire l’équation aux dérivées partielles de conduction
en un système différentiel ordinaire. Puis la technique des trapèzes généralisées est
utilisée pour les intégrations numériques dans le temps. Il ressort de l’analyse de
stabilité de l’algorithme de résolution du problème direct qu’il existe une relation
étroite entre le pas spatial de discrétisation et le pas d’intégration.
La défini-positivité de la matrice caractérisant le système linéaire du problème in-
verse assure qu’on peut identifier simultanément les deux paramètres thermophy-
siques.
Le champ de température et le vecteur flux de chaleur qui constituerait les donnés
d’entrée pour la méthode inverse sont générés dans les mêmes conditions que les
données d’entrée dans le cadre expérimental.
Les analyses statistiques des résultats pour divers pas spatiaux ont permis de choisir
un pas spatial proche d’un pixel. Puisque les valeurs des coefficients de conductivité
thermique des polymères sont faibles, nous choisirons un pas de temps supérieur
ou égale au pas de temps lors de la simulation numérique. Le fait que les valeurs
d’entrée pour la simulation se trouvent dans les intervalles de confiance permet de
valider l’algorithme.
Le coût de calcul de cet algorithme est très réduit par rapport à l’algorithme du
gradient conjugué.
Dans l’algorithme présenté pour identifier les paramètres thermophysiques, les vec-
teurs température sont directement incorporés dans le système différentiel ordinaire
obtenu par la MEF
V.8. CONCLUSION 119
L’algorithme d’identification est basé sur des techniques des moindres carrés.
Notre méthode des moindres carrés n’as pas consisté à minimiser des carrés des
différences des solutions obtenues par séries et des solutions obtenues expérimenta-
lement. Nous avions cherché directement à déterminer la conductivité et la chaleur
volumique qui rendent le système linéaire résoluble. Faisant ainsi, nous avons utilisé
la linéarité de ces paramètres. L’équation de conduction de la chaleur et sa forme
discrétisée par éléments finis sont linéaires par rapport à la conductivité thermique
et la chaleur volumique. Cet algorithme d’identification sera exploité dans le cha-
pitre suivant pour identifier les paramètres thermophysiques k et ρc de polymères à
partir de mesure de champs de température expérimentaux.
120 V. Construction d’un Algorithme d’Identification de Paramètres Thermophysiques
Chapitre VI
Identification de Paramètres
Thermophysiques de Polymères
Sommaire
VI.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 122
VI.2 Méthodes de mesures de champ de température . . . . . 123
VI.3 Dispositifs expérimentaux . . . . . . . . . . . . . . . . . . 124
VI.3.1 Mesure de la température par caméra infrarouge . . . . . 124
VI.3.2 Conditions expérimentales . . . . . . . . . . . . . . . . . . 125
VI.4 Montage expérimental . . . . . . . . . . . . . . . . . . . . 126
VI.5 Polymères translucide, semi-amorphe et amorphe étudiés127
VI.6 Traitement des images infrarouges . . . . . . . . . . . . . 128
VI.6.1 Historique des types de lissage . . . . . . . . . . . . . . . 128
VI.6.2 Logiciel utilisé . . . . . . . . . . . . . . . . . . . . . . . . 129
VI.7 Exploitation des données expérimentales . . . . . . . . . 129
VI.7.1 Calcul du flux de chaleur linéique . . . . . . . . . . . . . . 130
VI.7.2 Paramètres thermophysiques du polyméthylméthacrylate . 131
VI.7.3 Paramètres thermophysiques du polyéthylène transparent 134
VI.7.4 Paramètres thermophysiques du polyéthylène amorphe . . 138
VI.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
122 VI. Identification de Paramètres Thermophysiques de Polymères
VI.1 Introduction
∂T ∂2T
ρc =k 2, x ∈ [0, L]. (VI.1)
∂t ∂x
A cette équation on adjoint les conditions aux limites de flux traduisant les
VI.2. MÉTHODES DE MESURES DE CHAMP DE TEMPÉRATURE 123
dθ
ρc + kKθ(t) = F.
dt
La caractérisation des polymères passent donc pour l’estimation des paramètres
thermophysiques k et ρc dans (V.7). Ceci revient à poser le problème inverse qui
rappelons-le, consiste à déterminer ρc et k à partir de champ de température mesuré
sur des polymères. Il a été montré dans le chapitre précédent que la détermination
de ces paramètres thermophysiques de l’équation (V.7) revient à la recherche du
vecteur u de (voir section V.6)
Au = f .
Nous utilisons une caméra infrarouge avec laquelle on peut avoir accès aux
champs de température locaux ou une cartographie complète du matériau.
La mesure par rayonnement infrarouge émis par un corps opaque placé dans un
environnement correspond :
– au rayonnement de l’objet qui est affecté par la traversée de l’atmosphère,
– au rayonnement émis par l’atmosphère elle-même,
– au rayonnement émis par l’environnement, réfléchi par l’échantillon et modifié
par la traversée de l’atmosphère.
VI.3. DISPOSITIFS EXPÉRIMENTAUX 125
Une caméra de type FLIR SC3000 à matrice de détecteurs QWIP refroidie par
cycle Stirling dont le domaine spectral est de 8 à 9 µm ([Link].2) est utilisées pour la
prise des images. D’après les données techniques, données vérifiées par les travaux en
température est mesurée sur les deux échantillons. La différence des mesures sur
les deux éprouvettes donne la température relative en procédant de telle façon que
l’éprouvette témoin ne soit soumise à aucun flux de chaleur pouvant faire varier sa
température.
Afin d’avoir les meilleurs résultats possibles, nous devons nous placer dans de
bonnes conditions expérimentales :
– l’échantillon est peint en noir mat afin de le considérer comme un corps noir
avec une émissivité voisine de 1.
– la caméra infrarouge est placée le plus proche possible de l’échantillon et per-
pendiculairement à celui-ci.
– des obstacles sont mis en place de part et d’autre de la caméra pour minimiser
les réflexions.
126 VI. Identification de Paramètres Thermophysiques de Polymères
Le bruit des images infrarouges (±0, 02◦ C sur la température relative) peut pro-
venir du type de détecteur, des réflexions parasites à la surface du matériau ou du
matériau lui-même. Le fait de calculer les dérivées successives engendre du bruit
supplémentaire. Un lissage devient donc indispensable. Il existe plusieurs méthodes
permettant d’améliorer les images afin d’obtenir les meilleures informations pos-
sibles.
Notre logiciel d’acquisition repose sur deux types de filtrage : filtrage passe-bas
et filtrage temporel. Ce logiciel de filtrage est développé au sein de l’équipe Pho-
tomécanique et Rhéologie du Laboratoire de Mécanique des Solides de l’Université
de Poitiers. Il permet d’avoir la position, les valeurs de la température, les dérivées
temporelles et spatiales (les gradients de température), le temps, etc.
Pour éliminer le bruit provenant des dérivées spatiales, un filtre passe-bas clas-
sique peut être appliqué sur la température et ses dérivées spatiales. En dimension
un, le principe du filtrage passe-bas de dimension ν est :
Pν
θ(i, j, t)
θi (t) = l=1 , j = 1, ..., ν ; i = 0, ..., n
ν
L’utilisation d’un filtrage passe-bas avec une taille de filtre trop grande entraîne un
écrasement des valeurs maximales.
P Ri2
qs = =
s πrN l
L’éprouvette constituée de deux plaques identiques est chauffée avec des intensi-
tés de courant différentes. La carte thermique (VI.6) prise au cours du temps donne
la température à la surface de l’éprouvette. La température initiale du spécimen
enregistrée par la caméra infrarouge est de l’ordre de 20◦ C.
132 VI. Identification de Paramètres Thermophysiques de Polymères
Fig. VI.8 – Profils de température dans le Plexiglass pour q = 256W/m ; seuls les champs de température relevés sur
la partie [0, 20] mm ont été exploités.
EXPLOITATION DES DONNÉES EXPÉRIMENTALES
133
134 VI. Identification de Paramètres Thermophysiques de Polymères
Comme le montre le tableau VI.1, les valeurs estimées pour les paramètres thermo-
physiques ρc et k du plexiglass sont concordantes avec les données du fabricant. La
sensibilité de ρc et k n’est pas importante à la variation de température. Ceci peut
confirmer l’hypothèse selon laquelle les paramètres thermophysiques ρc et k pour de
faibles variations de température sont des constantes.
Les résultats obtenus sur trois essais sur polyéthylène transparent sont dans le ta-
bleau VI.2.
Pour une quantité de chaleur linéique q = 104, 15 W/m, les valeurs estimées des pa-
ramètres thermophysiques sont très mauvaises. Nous expliquons cette situation par
la très faible variation de la température qui rend le système (V.17) non résoluble.
En effet, d’une part, la première composante A11 de la matrice A est calculée avec
le vecteur gradient de température :
Z tf
dθ dθ
A11 = hC , C idt.
0 dt dt
Le profil du gradient de température est donné par la figure VI.11 ; le produit scalaire
hCdθ/dt, Cdθ/dti est de l’ordre de 10−6 , donc la composante A11 presque nulle.
D’autre part, nous avons
Z tf
A22 = hKθ(t), Kθ(t)idt,
0
Fig. VI.11 – Profils du gradient de température du PET transparent pour q = 104, 83W/m ; seuls les champs de
température relevés sur la partie [0, 20] mm ont été exploités.
EXPLOITATION DES DONNÉES EXPÉRIMENTALES
137
138 VI. Identification de Paramètres Thermophysiques de Polymères
Des tableaux VI.1, VI.2 et VI.3, nous pouvons résumer les incertitudes relatives
obtenues sur les paramètres identifiés dans le tableau VI.4
PMMA ∆k/k 2, 7% 1, 2%
∆ρc/ρc 3, 01% 2, 6%
PET transparent ∆k/k 1, 8% 2, 9% 2, 3%
∆ρc/rhoc 2% 19, 5% 3, 2%
PET Amorphe ∆k/k 2, 8% 3% 1, 6%
∆ρc/rhoc 3, 3%% 3, 2% 1, 5%
l’ordre de 3%. Les incertitudes relatives calculées sur k sont beaucoup plus faibles
dans le cas expérimental que dans le cas simulé. Selon une loi édictée par Fourier,
le flux de chaleur est proportionnel au gradient de température et le coefficient
de proportionnalité est égal à la conductivité thermique. Le bruit extrême sur les
données de température non réalistes dans la simulation pourrait se répercuter sur
la conductivité therùique et il en est de même pour le buit sur le flux de chaleur.
VI.8 Conclusion
Dans nos travaux sur les problèmes inverses en mécanique et en thermique, nous
avons présenté les problèmes d’identification de paramètres. Nous avons traité en
particulier l’identification des caractéristiques d’inertie d’un solide rigide en méca-
nique et l’identification de la conductivité thermique et de la chaleur volumique de
matériaux homogènes isotropes en thermique.
Dans le cas de la mécanique, les caractéristiques d’inertie du solide interviennent
linéairement dans les équations du mouvement sans aucune approximation. En ther-
mique, la chaleur volumique intervient linéairement ainsi que la conductivité ther-
mique dans le cas de la loi linéaire isotrope de Fourier.
La procédure qui a conduit aux identifications se résume en quatre étapes :
(i) Etudier théoriquement le caractère identifiable des paramètres mécaniques ou
thermophysiques du problème inverse.
(ii) Tenir compte des propriétés des matrices et vecteurs décrivant le problème
inverse : si les paramètres à identifier se regroupent dans une matrice symétrique
et positive, il faut tenir compte de la positivité de la matrice dans le processus de
minimisation par la technique des moindres carrés.
(iii) Valider le modèle inverse à partir des simulations des données conformes à
des conditions expérimentales.
(iv) Filtrer les mesures expérimentales afin de ne retenir que des matrices et des
vecteurs de mesure peu bruités.
L’application de ces démarches nous a permis d’identifier les caractéristiques
d’inertie d’un solide rigide. L’identification de ces paramètres est facilitée par la
formulation matricielle des équations du mouvement du solide. Les caractéristiques
d’inertie regroupées dans une matrice 4 × 4 symétrique (tenseur d’inertie de Poinsot
généralisé) sont identifiées par application d’un algorithme du type gradient conjugué
projeté.
144
tion :
1. Modélisation du phénomène mécanique ou thermique
2. Etude et construction d’un algorithme de résolution du problème direct
3. Etude et construction d’un algorithme de résolution du problème inverse
4. Validation de l’algorithme du problème inverse à partir de simulations
5. Application de l’algorithme du problème inverse à des mesures expérimentales.
En mécanique les perspectives concernent l’extension d’identification des carac-
téristiques d’inertie d’un solide à l’identification des caractéristiques d’inertie d’un
système multicorps tel qu’un robot ou un athlète.
En thermique, toujours dans le cadre des lois de Fourier linéaires, nous envisageons
d’étendre notre approche d’identification de paramètres thermophysiques des ma-
tériaux homogènes et isotropes aux matériaux homogènes non isotropes, puis aux
matériaux non homogènes et non isotropes. Nous pensons que cette future approche
pourra se tranférer directement à la caractérisation des capteurs hybrides photovol-
taïques intégrés aux batiments.
146
BIBLIOGRAPHIE 147
Bibliographie
[9] Anis Ammous, Faycal Sellami, Kaicar Ammous, Herve Morel, Bruno Al-
lard et Jean-Pierre Chante : Developing an equivalent thermal model for
discrete semiconductor packages. International Journal of Thermal Sciences,
42(5):533–539, mai 2003.
[10] Dominique Baillis et Jean-Fran¸cois Sacadura : Thermal radiation proper-
ties of dispersed media : theoretical prediction and experimental characteriza-
tion. Journal of Quantitative Spectroscopy and Radiative Transfer, 67(5):327–
363, décembre 2000.
[11] D. Balageas, L. : Nouvelle méthode d’interprétation des thermogrammes
pour la détermination de la diffusivité thermique par la méthode impulsionnelle
(méthode "flash"). Revue [Link]., 17(4):227–237, 1982.
[12] D. L. Balageas : Nouvelle méthode d’interprétation des thermogrammes
pour la détermination de la diffusivité thermique par la méthode impulsionelle.
Rev. Phys. Appl., 17:227–237, 1982.
[13] A.G. Bell : Upon the production of sound by radiant energy. Phil Mag,
11:510 – 528, 1881.
[14] Rachid Bennacer, Abdelwahab Tobbal, Hassen Beji et Patrick Vasseur :
Double diffusive convection in a vertical enclosure filled with anisotropic porous
media. International Journal of Thermal Sciences, 40(1):30–41, janvier 2001.
[15] MaItine Bergounioux et Fulbert Mignot : Contrôle dans les inéquations
et multiplicateurs de lagrange. Comptes Rendus de l’Académie des Sciences -
Séries I - Mathematics, 329(7):607–612, octobre 1999.
[16] S. Boyd et L. Vandenberghe : Convex Optimization. Cambridge University
Press, 2004.
[17] Haïm Brezis : Analyse fonctionnelle : Théories et applications. 1983.
[18] Jean Cea : Optimisation : théories et algorithmes. DUNOD, 1971.
[19] F. Cernuschi, P.G. Bison, A. Figari, S. Marinetti et E. Grinzato :
Thermal diffusivity measurements by photothermal and thermographic tech-
niques. International Journal of Thermophysics, 25(2):439–457, 2004.
[20] Karine Champion-Réaud : Méthode d’étude des vibrations d’un système
mécanique non basée sur le calcul de ses modes propres. Thèse de doctorat,
École Normale Supérieure de l’Aéronautique et de l’Espace, 2002.
BIBLIOGRAPHIE 149
[21] H.-T. Chen, S.-Y. Lin, H.-R. Wang et L.-C. Fang : Estimation of two-sided
boundary conditions for two-dimensional inverse heat conduction problems.
International Journal of Heat and Mass Transfer, 45(1):15–23, 2002.
[22] Y. Chen et J.E. McInroy : Estimating symmetric, positive definite ma-
trices in robotic control. In Proceedings - IEEE International Conference on
Robotics and Automation, volume 4, pages 4269–4274, Dept. of Comp. Sci.
and Engineering, Penn. State University, University Park, PA 16801, United
States, 2002.
[23] André Chrysochoos* et Robert Peyroux : Analyse expérimentale et modé-
lisation numérique des couplages thermomécaniques dans les matériaux solides.
Rev. Gén. Therm., 37:582–606, 1998.
[24] King-wah Eric Chu : Symmetric solutions of linear matrix equations by matrix
decompositions. Linear Algebra and its Applications, 119:35–50, juillet 1989.
[25] P. G. Ciarlet : Introduction a l’analyse numérique matricielle et à l’optimi-
sation. MASSON, 1982.
[26] Damien Claire, Stephane Roux et Fran¸cois Hild : Identification de conduc-
tivités thermiques et de propriétés élastiques locales par analyse de champs :
Identification of thermal conductivities and elastic properties by field analyses.
Mecanique & Industries, 4(6):655–665, 2003.
[27] R. Coquard, D. Baillis et D. Quenard : Experimental and theoretical
study of the hot-wire method applied to low-density thermal insulators. Inter-
national Journal of Heat and Mass Transfer, 49(23-24):4511–4524, novembre
2006.
[28] S. André D. Maillet et A. Degiovanni : Les erreurs sur la diffusivité
thermique mesurée par méthode flash : confrontation théorie-expérience. J.
Phys. III France, 3:883–809, 1993.
[29] B. Daccorogna : Introduction au calcul des variations. Presses Polytech-
niques et Universitaires Romandes, 1992.
[30] A. J. Davies : Finite Element methods : A First Approch. Clarendon Press.
Oxford, 1980.
[31] A. Degiovanni : Correction de longueur d’impulsion pour la mesure de la
diffusivité thermique par méthode flash. International Journal of Heat and
Mass Transfer, 30(10):2199–2200, octobre 1987.
150 BIBLIOGRAPHIE
[43] Cristina Ferreira et Fernando C. Silva : Inertia theorems for pairs of ma-
trices. Linear Algebra and its Applications, 381:37–52, avril 2004.
[44] Cristina Ferreira et Fernando C. Silva : Inertia theorems for pairs of ma-
trices, ii. Linear Algebra and its Applications, 413(2-3):425–439, mars 2006.
[45] Ali Fguiri, Naouel Daouas, Naceur Borjini, M-Sassi Radhouani et Ha-
bib Ben Aissia : Experimental inverse analysis for the determination of boun-
dary conditions in the parallel hot wire technique. Experimental Thermal and
Fluid Science, 31(3):209–220, janvier 2007.
[54] Claire Gobbé, Sébastien Isernaa et Bruno Ladevie : Hot strip method :
application to thermal characterisation of orthotropic media. International
Journal of Thermal Sciences, 43(10):951–958, octobre 2004.
[55] G Golub, H. et F. Van Loan, C. : Matrix Computations, Third Edition.
The Johns Hopkins University Press, 1996.
[56] C. Goncalves, V., O. L. VILARINHO, A. SCOTTI et GUIMARÃES :
Estimation of heat ant thermal efficiency in gtaw process bys using inverse
techniques. Journal of Materials Processing Technology, 172:42–51, 2006.
[57] M.D. Greenberg : Foundation of applied mathematics. CPrentice-Hall, Inc.,
Englewood Cliffs, New York 07632, 1978.
[58] Allaire Gregoire : Analyse numérique et optimisation. 2005.
[59] E. Hensel : Inverse Theory and Applications for Engineers. Prentice Hall,
Englewood Cliffs, New Jersey 07632, 1990.
[60] Nicholas J. Higham : Computing a nearest symmetric positive semidefinite
matrix. Linear Algebra and its Applications, 103:103–118, mai 1988.
[61] Nicholas J. Higham : Symmetric procrustes problem. BIT (Copenhagen),
28(1):133–143, 1988.
[62] Dai Hua : Completing a symmetric 2 x 2 block matrix and its inverse. Linear
Algebra and its Applications, 235:235–245, mars 1996.
[63] Cheng-Hung Huang et Yan Jan-Yuan : An inverse problem in simulta-
neously measuring temperature-dependent thermal conductivity and heat ca-
pacity. International Journal of Heat and Mass Transfer, 38(18):3433–3441,
décembre 1995.
[64] Peter Ireman et Quoc-Son Nguyen : Using the gradients of temperature
and internal parameters in continuum thermodynamics. Comptes Rendus Me-
canique, 332(4):249–255, avril 2004.
[65] Kazufumi Ito et Karl Kunisch : Parabolic variational inequalities : The
lagrange multiplier approach. Journal de Mathematiques Pures et Appliques,
85(3):415–449, mars 2006.
[66] Y. Jarny, M. N. Ozisik et J. P. Bardon : A general optimization method
using adjoint equation for solving multidimensional inverse heat conduction.
International Journal of Heat and Mass Transfer, 34(11):2911–2919, novembre
1991.
BIBLIOGRAPHIE 153
[105] Yongxin Yuan et Hua Dai : Inverse problems for symmetric matrices with a
submatrix constraint. Applied Numerical Mathematics, –:–, –.
[106] S. Zagatti : On the minimum problem for nonconvex scalar functionals.
SIAM Journal on Mathematical Analysis, 37(3):982–995, 2006.
[107] Xing Zhang, Hua Gu et Motoo Fujii : Effective thermal conductivity and
thermal diffusivity of nanofluids containing spherical and cylindrical nanopar-
ticles. Experimental Thermal and Fluid Science, In Press, Corrected Proof:–
4524.
ANNEXE A Calcul des Composantes
des Matrices C et K Issues de la
Discrétisation par Eléments Finis de
l’Equation de Conduction de la
Chaleur
∂T (x, t) ∂ 2 T (x, t)
ρc =k (VI.5)
∂t ∂x2
avec condition initiale
T (x, 0) = T0 (x) donnée (VI.6)
1
ζ x + 1 − i si (i − 1)ζ ≤ x ≤ iζ
φi (x) = − ζ1 x + 1 + i si iζ ≤ x ≤ (i + 1)ζ
0
sinon
1 x − (n − 1) si (n − 1)ζ ≥ x ≥ nζ
φn (x) = ζ
0 sinon
Pour obtenir les matrices C et K et le veteur F, on multiplie directement l’équa-
tion aux dérivées partielles VI.5 par φi (x) puis on l’intègre sur [0, l].
Z l Z l 2
∂T (x, t) ∂ T (x, t)
ρc φi (x)dx = k φi (x)dx (VI.8)
0 ∂t 0 ∂x2
R l 2 T (x,t)
Calculons alors pour tout i compris entre 0 et n, l’intégrale k 0 ∂ ∂x 2 φi (x)dx
Le vecteur temérature θ dans le cadre unidimensionnel a pour composantes :
θi (t) = T (xi , t) θi (0) = T0 (xi ) i = 0, ..., n.
Pour i = 0, par intégration par partie on obtient :
Z l Z l 2
∂T (x, t) ∂ T (x, t)
ρc φ0 (x)dx = k φ0 (x)dx
0 ∂t 0 ∂x2
Z ζ 2
∂ T (x, t)
=k φ0 (x)dx
0 ∂x2
∂T (0, t) k
= −k + (θ1 (t) − θ0 (t))
∂x ζ
De même pour 1 ≤ i ≤ n − 1, l’intégration par parties conduit à :
Z l Z l 2
∂T (x, t) ∂ T (x, t)
ρc φi (x)dx = k φi (x)dx
0 ∂t 0 ∂x2
Z (i+1)ζ 2
∂ T (x, t)
=k φi (x)dx
(i−1)ζ ∂x2
k k k
= θi−1 (t) − 2 θi (t) + θi−1 (t)
ζ ζ ζ
Enfin pour i = n, on obtient :
Z l Z l 2
∂T (x, t) ∂ T (x, t)
ρc φn (x)dx = k φn (x)dx
0 ∂t 0 ∂x2
Z l
∂ 2 T (x, t)
=k φn (x)dx
(n−1)ζ ∂x2
∂T (l, t) k
=k + (θn−1 (t) − θn (t))
∂x p
ANNEXE 159
Z l Z l
T (x, t)φ0 (x)dx = T (x, t)φ0 (x)dx
0 0
Z ζ Z ζ
2
= θ0 (t) (φ0 ) (x)dx + θ1 (t) φ1 (x)φ0 (x)dx
0 0
ζ
= [2θ0 (t) + θ1 (t)]
6
Z l Z (i+1)ζ
T (x, t)φi (x)dx = T (x, t)φi (x)dx
0 (i−1)ζ
ζ
= [θi−1 (t) + 4θi (t) + θi+1 (t)]
6
160
Z l Z l
T (x, t)φn (x)dx = T (x, t)φn (x)dx
0 (n−1)ζ
ζ
= [2θn (t) + θn−1 (t)]
6
Rζ
0
T (x, t)φ0 (x)dx
R 2ζ
T (x, t)φ1 (x)dx
0
..
.
..
∂ R .
(i+1)ζ
∂t (i−1)ζ T (x, t)φi (x)dx
..
.
..
.
Rl
(n−1)ζ
T (x, t)φn (x)dx
2 1 0 0 ··· ··· ··· ··· 0 θ0 (t)
1 4 1 0 ··· ··· ··· ··· 0 θ1 (t)
. ... ... ... ... ... ... . . ..
0 .. . .
.
. ... ... .. .. .. .. . . .. ..
. . . . . . . .
ζ d ... . . . . . . .. .. .. .. . . .. ..
= ρc . . . . . . .
6 dt
. .
.. . . . . . . . . . . . . . . . . . . . . .. ..
.
.
.. . . ..
. ... ... ... .. .. ..
.
. . . 0
.
0 ··· ··· ··· ··· 0 1 4 1 θn−1 (t)
0 ··· ··· ··· ··· 0 0 1 2 θn (t)
ANNEXE 161
Notons par δθ[1 ... 1]T l’incertitude sur le vecteur de température θ(t) et par
δ θ̇[1 ... 1]T l’incertitude sur le vecteur dérivé de température θ̇(t) à tout instant t.
Calcul des incertitudes sur la matrice A et le vecteur f .
Sur la composante A11 , x = y = C θ̇(t) et xj = yj = C θ̇(tj ), j = 1, ..., m − 1.
j j
δhx , x i = δ hC θ̇(tj ), C θ̇(tj )i
1
.
C . , θ̇(tj )i
= 2δ θ̇hC
.
1
Pour les composantes A12 et A21 qui sont égales, posons xj = C θ̇(tj ) et yj = Kθ(tj )
δhxj , yj i = δ hC θ̇(tj ), Kθ(tj )i
1 1 0
. . .
= δ θ̇hK C . , θ(tj )i car K . =
.
. . .
1 1 0
Les composantes des matrices C et K ont été calculées dans la section V.3.1. On
164
√ um−1
v
δ θ̇ 2 u X
δA12 = δA21 = t (hw, θ(tj )i)2 (tj+1 − tj )2
tf 2 j=1
δA22 = 0 puisque le noyau de K contient les vecteurs colinéaires au vecteur [1 1 ... 1]T .
Calculons les incertitudes sur les composantes du vecteur f .
Rt
De l’expression f1 = 0 f hC θ̇(t), [q 0 ... 0]T idt, posons xj = C θ̇(tj ) et yj = [q 0 ... 0]T
Rt
et de de f2 = L2 /tf 0 f hKθ, [q 0 ... 0]T i, posons xj = Kθ(tj ) et yj = [q 0 ... 0]T .
En procédant de la même manière que pour les composantes de la matrice A, nous
obtenons
√ um−1 h
v
2u X i
δf1 = t qδ θ̇(hCr, e1 i)2 + δq(hθ̇(tj ), Ce1 i)2 (tj+1 − tj )2
2 j=1
√ 2 um−1
v
2 L δq u X
δf2 = t (hθ(tj ), Ke1 i)2 (tj+1 − tj )2
2 tf j=1
1 1
0 1
. 1
e1 = et r =
.
.
. .
0 1