La modélisation unidimensionnelle

La thermique

Comme nous l’avons montré dans le paragraphe (2.3) page 35 , les échanges thermiques au cours du fibrage sont de différentes natures. Compte tenu du caractère semi-transparent du verre, les échanges par rayonnement sont plus complexes à modéliser que s’il s’agissait d’un corps noir. 

 Influence des différentes contributions (convection et rayonnement) sur le profil de température

 Nous allons ici chercher à mettre en évidence la nature des échanges thermiques au cours du fibrage. Pour cela, nous allons repartir des résultats obtenus avec le modèle stationnaire et résoudre une équation thermique avec des contributions différentes. Dans le paragraphe précédent, nous sommes partis du profil de température expérimental afin de déterminer la géométrie de la fibre le long du chemin de fibrage. Nous allons à présent partir du profil géométrique de la fibre et résoudre l’équation de la chaleur pour retrouver le profil de température. L’équation thermique que nous allons résoudre est de la forme :Afin de tester l’influence de contributions comme la convection forcée ou le rayonnement, nous allons exprimer hT comme la somme d’un terme de convection forcée hc et d’un terme de rayonnement hr et résoudre l’équation (4.7) dans trois configurations :  Terme de convection forcée uniquement  Terme de rayonnement uniquement  Somme des deux termes La figure (4.4) montre que le profil de température le plus proche du profil réel est celui qui a été calculé avec une contribution de convection seule. Il semble que la contribution d’un terme de rayonnement (tel qu’il a été exprimé pour cette étude) fasse s’éloigner le profil calculé du profil réel.

Etude de stabilité linéaire 

Définition

 Il s’agit d’étudier l’évolution d’une petite perturbation de la solution stationnaire des équations du modèle. Typiquement, la forme perturbée d’une solution stationnaire f0(z) est décrite par l’expression : f(z, t) = f0(z) + ˜f(z) e λt. (λ complexe et ˜f à valeurs complexes). Cette forme perturbée est introduite dans les équations instationnaires du modèle qui sont ensuite linéarisées en ne gardant que les termes du premier ordre. Nous obtenons alors une relation entre λ et les paramètres du procédé, ce qui permet de calculer λ qui est appelée valeur propre du système. C’est la valeur propre de plus grande partie réelle qui permet de conclure sur la stabilité. Deux cas se présentent alors : si la partie réelle de cette valeur propre est négative alors le système tend vers sa solution stationnaire, nous sommes dans un cas stable. Si au contraire, la partie réelle est positive, la solution diverge de sa solution stationnaire et nous sommes alors dans un cas instable. Nous allons chercher la solution des équations instationnaires sous la forme : u(z, t) = u1(z) + ˜u(z) e λt (4.11) s(z, t) = s1(z) + ˜s(z) e λt (4.12) F(z, t) = F1(z) + F˜(z) e λt (4.13) T(z, t) = T1(z) + T˜(z) e λt (4.14) O`u u1, s1, F1 et T1 sont les solutions du système stationnaire (équations 4.1 à 4.4). Les fonctions ˜u, ˜s, F˜ et T˜ sont les fonctions perturbées inconnues également appelées fonctions propres du système. 

 La résolution

 La résolution de ce système a été faîte par deux méthodes qui vont s’avérer complémentaires l’une de l’autre. La première, dite méthode matricielle, va nous permettre d’obtenir un certain nombre de valeurs propres mais de manière assez peu précise. La seconde méthode, dite méthode de résolution de l’équation aux valeurs propres, est plus précise mais nécessite d’avoir une petite idée de la valeur propre cherchée.

Résolution par la méthode matricielle

Comme son nom l’indique, nous allons transformer le système afin de le mettre sous forme matricielle AX = λBX. Pour cela, nous utilisons les approximations des différences finies pour remplacer les dérivées et transformer ainsi le système différentiel en système algébrique. Nous discrétisons donc l’espace z ∈ [0;L] en N intervalles [zi ; zi+1]i=0..(n−1) et sur chacun de ces intervalles, nous approchons les dérivées par :    df dz ‘ 3f(z) − 4f(z − h) + f(z − 2h) h h = zi+1 − zi i = 1..(n − 1) Ce schéma n’étant utilisable qu’à partir du second point de discrétisation, nous utiliserons pour le premier point la forme plus classique : df dz ‘ f(z) − f(z − h) h Compte tenu des conditions aux limites, nous pouvons remarquer que : X s˜i et T˜ i sont définies pour i ∈ [1; n] car les conditions aux limites en z = 0 sont nulles. 61 Chapitre 4. La modélisation unidimensionnelle X F˜ i est définie pour i ∈ [0; n] car F˜ 0 = F1(0) non nulle. X u˜i est définie pour i ∈ [1; n − 1] car les valeurs en z = 0 et z = L sont nulles. Nous obtenons alors un système de 4n équations à 4n inconnues et fonction du paramètre λ, inconnu lui aussi. La résolution de ce système se fait numériquement, grâce à un sous-programme de la librairie IMSL utilisable en Fortran. Elle retourne comme résultats un vecteur contenant 4n valeurs propres complexes. Nous voyons donc clairement les limites de cette méthode. En effet, le nombre de valeurs propres obtenues est directement lié au nombre de points de discrétisation (puisque celui-ci implique la dimension de la matrice). Donc, il est clair que toutes les valeurs propres obtenues ne sont pas significatives. Il faut donc les trier. Cette méthode est d’autant plus mauvaise qu’elle est lourde à mettre en œuvre. En effet, le remplissage des matrices est fastidieux comme nous pouvons l’observer dans le détail joint en annexe A page 101. 

Résolution de l’équation aux valeurs propres

 Le système différentiel (équations (4.19) à (4.22)) est résolu par un algorithme de Runge et Kutta d’ordre deux en prenant en compte les conditions aux limites à une extrémité. Nous obtenons les fonctions propres solutions pour une valeur de λ donnée. Il reste à satisfaire la condition à l’autre bout : ˜u(L) = 0. Nous obtenons ainsi une équation aux valeurs propres de la forme f(λ) = 0 qui est résolue par une méthode de Newton. Il existe une infinité de valeurs de λ solutions de cette équation, cependant, la valeur propre à considérer sera celle de plus grande partie réelle car elle donne la stabilité. Il y a deux moyens d’aborder cette résolution, en considérant soit les conditions limites en haut du chemin de fibrage (en z = 0) soit en bas (en z = L).