Ecole Mohammadia d’Ingénieurs- GPI 2

Partie I : Initiation-Matlab
Quelques applications pour la filière Génie des Procédés
Industriels

Sommaire
1- Initiation

2- Création des fonctions

3- Représentation graphique des fonctions (2D et 3D)

4- Quelques Méthodes de Résolution numérique des équations non-linéaires

5- Quelques méthodes d’optimisation des systèmes complexes

6- Exemple d’application :

 Calcul des réacteurs chimiques

 Modélisation et simulation d’un collecteur solaire cylindro-parabolique
 Dimensionnement d’un évaporateur à film tombant

1
 Ecole Mohammadia d’Ingénieurs- GPI 2

1- Initiation
Le logiciel MATLAB (langage développé par la société The MathWorks) est un langage
qui intègre à la fois le calcul, la visualisation graphiques et la programmation dans un
environnement facile à utiliser. La structure de données de base est le tableau (calcul matriciel),
cette structure ne nécessite pas de dimensionnement cela permet d’utiliser des notations
mathématiques simples et familières et de résoudre beaucoup de problèmes techniques comme la
résolution numériques des équations non-linéaires et différentielles, traitement des problèmes
complexes d’optimisation …etc. Aussi, il fournit de nombreuses fonctions préprogrammées
regroupées en boîtes à outils (toolbox) pour de nombreux domaines d’application (comme :
signal processing, statistics, control theory, optimization, ...).

Démarrage de Matlab :

Lorsque vous lancez MATLAB pour la première fois, l’écran ressemble à celui de la
Figure ci-dessous :

2
 Ecole Mohammadia d’Ingénieurs- GPI 2

Opérations de calcul :

sin (Sinus) ⇒ >> sin(7)

+ (addition) ⇒ >> 7+3
cos (cosinus) ⇒ >> cos(7)
- (soustraction) ⇒ >> 7-3
log (log népérien) ⇒ >> log(7)
* (multiplication) ⇒ >> 7*3
log10 (log décimal ) ⇒ >> log10(7)
/ (division) ⇒ >> 7/3
tan (tangente) ⇒ >> tan(7)
^ (puissance) ⇒ >> 7^3
exp (exponentielle) ⇒ >> exp(7)

Affectation des variables

On peut spécifier le nom d’une variable pour stocker le résultat. Et on pourra l’utiliser
plus tard dans d’autres opérations comme la montre la figure en haut. L’ordre de priorité des
opérations est le suivant :

1 : les parenthèses ( ), 2 : la puissance ^, 3 : multiplication et division * / et 4 : somme et

soustraction + et –

Exemples :

3
 Ecole Mohammadia d’Ingénieurs- GPI 2

Construction des matrices :

On sépare les éléments du vecteur par une virgule ou par un espace pour obtenir un
agencement en ligne, exemples :

On sépare les éléments par un point-virgule pour obtenir un agencement en colonne

4
 Ecole Mohammadia d’Ingénieurs- GPI 2

Pour construire une matrice : les éléments d’une même ligne sont séparés par des virgules
et les lignes de la matrice sont séparées par un point-virgule

5
 Ecole Mohammadia d’Ingénieurs- GPI 2

Opérations et calcul matriciel

Lorsque on utilise un ‘ .’ avant les opérations de multiplication, division et puissance, le calcul

s’effectue par élément à l’intérieur de la matrice (voir les exemples : a*d et a.*d, d=a^-1 et d=a.^-1
…etc)

6
Ecole Mohammadia d’Ingénieurs- GPI 2

7
 Ecole Mohammadia d’Ingénieurs- GPI 2

2. Création des fonctions

Nous avons mentionné précédemment qu’il y a des fonctions MATLAB préprogrammées comme
les fonctions max(), max(max()), min(), size(), clear … etc. Matlab donne la possibilité de créer ses
propres fonctions au sein d’un M-file en utilisant des scripts. Ces fonctions peuvent accepter des données
en entrée et renvoyer des données en sortie. Chaque fonction a son propre espace en mémoire et ne voie
pas les données du workspace. Seules les données qu’on aura spécifiées comme arguments d’entrée
seront utilisables par la fonction.

Ouverture d’un nouveau script

8
 Ecole Mohammadia d’Ingénieurs- GPI 2

Exemple : Création de la fonction nommée my_function(x)

( ) si x=<0

( ) si x>0

9
 Ecole Mohammadia d’Ingénieurs- GPI 2

3-Représentation graphique

10
Ecole Mohammadia d’Ingénieurs- GPI 2

11
Ecole Mohammadia d’Ingénieurs- GPI 2

12
Ecole Mohammadia d’Ingénieurs- GPI 2

13
 Ecole Mohammadia d’Ingénieurs- GPI 2

4- Quelques Méthodes de Résolution numérique des équations non-linéaires

Exemple 1 :

% utilisation de la méthode de Newton-Raphson

% la fonction à résoudre f(x)=5*x^3 + 3*exp(6*x-2)-5
% la fonction dérivée de f(x)est f'(x)=15*x^2+18*exp(6*x-2)
% le calcul itératif se base sur la relation: Xn+1=Xn - f(Xn)/f'(Xn)

% Initialisation de la solution X0

X0=1;
Error=1000;
% boucle While
while Error>0.0000000001
x=X0;
F_x=5*x^3 + 3*exp(6*x-2)-5;
dF_x=15*x^2+18*exp(6*x-2);
Xn=x-(F_x/dF_x);
Error=abs(Xn-x);
X0=Xn;
end

% boucle for
X1=1;
for i=1:3000
x1=X1;
F_x1=5*x1^3 + 3*exp(6*x1-2)-5;
dF_x1=15*x1^2+18*exp(6*x1-2);
X1n=x1-(F_x1/dF_x1);
X1=X1n;
end

% Verification graphique
y=0.404:0.0001:0.41;
F_y=5*y.^3 + 3*exp(6*y-2)-5;
plot(y,F_y,'--x')
grid
xlabel('Variable X')
ylabel('F(X)')

Résultat :

14
 Ecole Mohammadia d’Ingénieurs- GPI 2

Exemple 2 :

% Méthode de Newton-Raphson
% Système de deux équations non linéaires
% F1=x^2 + y - 1
% F2=x^3 - y
% DF1_x=2*x, DF1_y=1, DF2_x=3*x^2, DF2_y=-1
% Matrice Jacobienne : J=[DF1_x , DF1_y ; DF2_x , DF2_y ]
% Solution Xn+1=[xn+1;yn+1]=[xn;yn]-J(xn,yn)^(-1)*[F1(xn,yn);F1(xn,yn)];

X0=[0.5;0.5];
Error=1000;
% boucle While
while Error>10^(-20)
X=X0;
F1=X(1,1)^2 + X(2,1)-1;
F2=X(1,1)^3 -X(2,1);
DF1_x=2*X(1,1);
DF1_y=1;
DF2_x=3*X(1,1)^2;
DF2_y=-1;
J=[DF1_x , DF1_y ; DF2_x , DF2_y ];
Xn=X-(J^-1)*[F1;F2];
Error=abs(Xn(1,1)-X(1,1))+1000*abs(Xn(2,1)-X(2,1));
X0=Xn;
end

Résultat et vérification de la solution

15
 Ecole Mohammadia d’Ingénieurs- GPI 2

Exemple 3 :

% Utilisation du Solveur Matlab

% Le système à résoudre est F1(X1, X2, X3), F2(X1, X2, X3) et F3(X1, X2, X3)

syms X1 X2 X3

F1=-exp(X1^5 -X2^2)+ X3^3-5;

F2=X1^3-X2- exp(X3+2);
F3=X1^2 + X3^2-4;

eqns = [F1,F2,F3];
S = solve(eqns,X1,X2,X3);
X1=S.X1;
X2=S.X2;
X3=S.X3;

Exemple 4 :

% L'équation à résoudre est: dy/dx=f(x,y)=f=3*X+Y^0.3-1;

% Condition initiale x0=0 , y0=2
% Trouver Y pour x=1, x=2, et x=3

% Méthode de Runge-Kutta : Equations différentielles

% Xn+1=Xn + h (h: Size Step)
% Yn+1=Yn+ h*T(Xn,Yn,h)

% T(Xn,Yn,h)=(1/6)*(k1+2*k2+2*k3+k4)
% k1=f(Xn,Yn)
% k2=f(Xn+h/2,Yn+(h/2)*k1)
% k3=f(Xn+h/2,Yn+(h/2)*k2)

16
 Ecole Mohammadia d’Ingénieurs- GPI 2

% k4=f(Xn+h,Yn+h*k3)

% Etape 1: Creation de la fonction f(x,y)= 3*X+Y^0.3-1 dans un script;

function f=f(X,Y)
f=3*X+Y^0.3-1;
end

% Etape 2: dans un autre script on applique la méthode de Runge-Kutta

h=0.01;
X=0:h:3;

dim=size(X);
n=dim(1,2);
X(1)=0;
Y=zeros(1,n);
Y(1)=2;
for i=2:n
k1=f(X(i-1),Y(i-1));
k2=f(X(i-1)+h/2,Y(i-1)+(h/2)*k1);
k3=f(X(i-1)+h/2,Y(i-1)+(h/2)*k2);
k4=f(X(i-1)+h,Y(i-1)+h*k3);
T=(1/6)*(k1+2*k2+2*k3+k4);
X(i)=X(i-1)+h;
Y(i)=Y(i-1)+h*T;
end

plot(X,Y,'-o')
xlabel('Variable X')
ylabel('Fonction Y')
grid

17
 Ecole Mohammadia d’Ingénieurs- GPI 2

5- Quelque méthode d’optimisation des systèmes complexes

Méthode du simplex : Problèmes linéaires simples avec des contraintes.

Fonction (préprogrammées) à utiliser : linprog(f,A,b,Aeq,beq,lb,ub)

Exemple 1:

La fonction à minimiser est : f(x,y)=3*x+6*y

Les contraintes d’inégalité sont définies comme suit :

Code :

% La fonction à minimiser est : f(x,y)=3*x+6*y

% Les contraintes d’inégalité sont définies comme suit :
% x+y=<0
% x-y=<5
% -4x-3y=< 3

f=[3 6];
A=[1 1; 1 -1; -4 -3];
b=[0 5 3];
x = linprog(f,A,b);

Résultat :

18
 Ecole Mohammadia d’Ingénieurs- GPI 2

Exemple 2:

% La fonction à minimiser est : f(x,y,z,w)=x+6*y -8*z

% Les contraintes d’inégalité sont définies comme suit :

% -5⫹x⫹5 et 0⫹y⫹15 et 0⫹z⫹15

% x+y-2*z⫹-10
% x+y+z⫹25
% 4x-3y⫹4

% Les contraintes d’égalité :

% x+2*y-z=15

f=[1 6 -8];
A=[1 1 -2; 1 1 1; 4 -3 0];
b=[-10 25 4];
Aeq=[1 2 -1];
beq=15;
lb=[-5 0 0];
ub=[5 15 15];
Solution= linprog(f,A,b,Aeq,beq,lb,ub);

Résultat:

19
 Ecole Mohammadia d’Ingénieurs- GPI 2

Méthode de l’algorithme génétique (Genetic Algorithm optimization method)

La méthode de l’algorithme génétique est une méthode qui permet de générer une solution rapide et
faisable pour des problèmes d’optimisation complexe et hautement non-linéaires.

Exemple :

Fonction à minimiser : ( ) ( )
( )
Contraintes d’inégalité :
( )
( ) ( )

Contraintes d’égalité :

La fonction à utiliser : ga(ObjectiveFunction,nvars,[],[],[],[],LB,UB,ConstraintFunction,options);

L’utilisation de cette fonction préprogrammée dans le toolbox de matlab nécessite d’abord la création de
deux fonctions dans deux scripts séparés : la fonction objective et la fonction des contraintes. Et
finalement, la création d’un troisième script pour l’exécution de la fonction ‘ga’.

La fonction objective :

function [f] = Fitness(x)

x1=x(1);
y=x(2);
z=x(3);
w=x(4);
t=x(5);
v=x(6);
f=exp(x1^2*y*t)+w/log(z)-w^3+z^0.5-t*v+v^4;
end

La fonction des contraintes :

function [c, ceq] = Constraint(x)

x1=x(1);
y=x(2);
z=x(3);
w=x(4);
t=x(5);
v=x(6);
c =[x1^3+y^0.5-log(w*t)-3;cos(x1)+sin(v)-0.9]; % contraintes d'Inégalité
ceq =x1+z^5+t^2-100; % contraintes d'égalité

20
 Ecole Mohammadia d’Ingénieurs- GPI 2

end

Le script d’exécution:

clear all
clc
nvars =6; % Number of variables
LB =ones(1,6); % Lower bound
UB =100*ones(1,6); % Upper bound
type Fitness;
X0 =[2,50,0,3,20,7] ; % Initialization
ObjectiveFunction = @(x) Fitness(x);
ConstraintFunction = @(x)Constraint(x);
[Link] = X0;

% Next we run the GA solver.

options = gaoptimset(options,'PlotFcn',{@gaplotbestf,@gaplotmaxconstr}, ...
'Display','iter','Useparallel','Always');
[x,fval] = ga(ObjectiveFunction,nvars,[],[],[],[],LB,UB, ...
ConstraintFunction,options);

Résultat :

21
 Ecole Mohammadia d’Ingénieurs- GPI 2

Autres méthodes que vous pouvez trouver sur le site [Link] ( Méthode de
Lagrange, PSO : Particle Swarm Optimization, Simulated annealing (SA), Neural Network … etc )

22
 Ecole Mohammadia d’Ingénieurs- GPI 2

6- Exemples d’application :
 Calcul des réacteurs (Exercices relatives au cours des réacteurs homogènes idéaux)

A-1- Exercice : Annexe-Page 10

Données :

 Réaction chimique : déshydrogénation en phase gazeuse du benzène :

( ) ( ) ( ) ; Réaction1 : ( )

( ) ( ) ( ) ( ) ; Réaction2 : ( )

B : Benzène, D : Diphényle, T : Tri-phényle, H : Hydrogène, p est la pression partielle en atm et

les constantes des vitesses de réaction sont données par (Température est exprimée en K) :

( ) et ( )

r1 est exprimée en lbmol de B réagit/( hr.ft3) et r2 est donnée en lbmol de T ou D réagit/(hr.ft3).

 Type du réacteur : Réacteur Tubulaire

Déterminer :

 La conversion totale du benzène (B) en D et T ;

 Le volume du réacteur nécessaire pour traiter 10000 lb/hr de B en fonction de la
conversion totale de B. (l’alimentation est constituée B pur) ;
 La vitesse spatiale qui permet d’obtenir une conversion maximale.

Solution : (hypothèses : réacteur piston, isotherme, régime permanent, gaz parfaits … etc.)

On définit le taux de conversion de B en D par :

( )

Et le taux de conversion de B en T :

( )

Equation bilan de matière par rapport à D (considérant un volume différentiel dV=[Link]):

23
 Ecole Mohammadia d’Ingénieurs- GPI 2

( ) => ( )
( ) ( )
Equation bilan de matière par rapport à T :

=> ( )
Expression des pressions partielles en fonction des conversions X1 et X2 :

( ) ( ) ; ( ) ( )

( ) ; ( ) ( )
Les vitesses des deux réactions sont exprimées par :

(( ) ( )( ) )

(( )( ) ( ) )
De ce fait, les équations des bilans de matière à résoudre sont exprimées par les deux équations
différentielles ci-dessous :

( ) ( ( ) ( )) ( )
( )

( ) ( ) ( )
( )

Résolution Numérique: Méthode de Runge-Kutta :

Conditions initiales : V/FBf=0, X1,0=0 et X2,0=0. On choisit un incrément de ∆(V/FBf)=0.01

Pour chaque itération n+1 :

∆( ) ( ( ))

∆( ) ( ( ))

∆( ) ( ( ))

∆( ) ( ( ))

∆( ) ( ( ))

∆( ) ( ( ))

24
 Ecole Mohammadia d’Ingénieurs- GPI 2

∆( ) ( ( ))

∆( ) ( ( ))

( ) ( )

( ) ( )

( ) ( ) ∆( )

Programmes à réaliser sur Matlab :

Création des deux fonctions f1(X1,X2) et f2(X1,X2) :

25
 Ecole Mohammadia d’Ingénieurs- GPI 2

Création du programme de résolution et affichage des résultats sous forme de figures:

T=1033.15; % Température en kelvin
k1=14.96*(10^6)*exp(-15200/T);
k2=8.67*(10^6)*exp(-15200/T);

size=40; % Nombre d'itérations

VF=zeros(1,size); % rapport entre volume du réacteur et débit d'alimentation
en B: V/FBf
X1=zeros(1,size);
X2=zeros(1,size);
DeltaVF=0.04; % Incrément pour V/FBf

for i=2:size
kk0=DeltaVF*f1(X1(i-1),X2(i-1),k1,k2);
m0=DeltaVF*f2(X1(i-1),X2(i-1),k2);

kk1=DeltaVF*f1(X1(i-1)+kk0/2,X2(i-1)+m0/2,k1,k2);
m1=DeltaVF*f2(X1(i-1)+kk0/2,X2(i-1)+m0/2,k2);

kk2=DeltaVF*f1(X1(i-1)+kk1/2,X2(i-1)+m1/2,k1,k2);
m2=DeltaVF*f2(X1(i-1)+kk1/2,X2(i-1)+m1/2,k2);

kk3=DeltaVF*f1(X1(i-1)+kk2,X2(i-1)+m2,k1,k2);
m3=DeltaVF*f2(X1(i-1)+kk2,X2(i-1)+m2,k2);

X1(i)=X1(i-1)+(1/6)*(kk0+(2*kk1)+(2*kk2)+kk3);
X2(i)=X2(i-1)+(1/6)*(m0+(2*m1)+(2*m2)+m3);
VF(i)=VF(i-1)+DeltaVF;

end

r1=k1.*(((1-X1-X2).^2)-(((X1/2-X2).*(X1/2+X2))./k1));
r2=k2.*((1-X1-X2).*(X1/2-X2)-((X2.*(X1/2+X2))./k2));

plot(VF,100*X1,'-x',VF,100*X2,'-+',VF,100*(X2+X1),'-*')
legend('X1','X2','Xtotal')
title('Taux de Conversion')
xlabel('rapport Volume/Débit molaire (ft^3/(lbmol/hr))')
ylabel('Conversion (%)')
grid

figure
plot(VF,r1,'-x',VF,r2,'-+')
title('Vitesse de réactions')
legend('r1','r2')
xlabel('rapport Volume/Débit molaire')
ylabel('Vitesse de réaction (lbmol/([Link]^3)')
grid

26
 Ecole Mohammadia d’Ingénieurs- GPI 2

Résultats obtenus :

A-2- Exercice : Annexe-Page 22

Données :

 Réaction chimique : décomposition d’huile castor en huile sèche:

() () ( ) ;

; Vitesse de réaction : avec C : concentration de l’acide dans l’huile

humide exprimée en kg d’acide/kg d’huile ; r exprimée en kg d’acide/ (s.m3).
La constante de la vitesse de réaction est donnée par (T en (K) et R en (cal/(mol.K))) :
( )

 Type du réacteur : réacteur à cuve idéal fermé et adiabatique

 Charge initiale : 227 kg d’huile de densité 0.9 à 340°C, contenant 0.156 kg d’acide par
kg d’huile.
 Chaleur de réaction (endothermique) (∆ )
 Capacité thermique moyenne du mélange réactionnel ( )

Travail demandé :

27
 Ecole Mohammadia d’Ingénieurs- GPI 2

Déterminer la variation du taux de conversion et de la température en fonction du temps de

réaction dans les deux cas suivants :

 Cas 1 : Réacteur adiabatique ;

 Cas 2 : Réacteur non adiabatique : chaleur fournie au réacteur : Q’=52700 J/s.
Représenter les résultats dans le même graphique.

Solution : (hypothèses: réacteur à cuve parfaitement agité, adiabatique, densité constante

du mélange réactionnel … etc.)

Relation entre concentration d’acide et taux de conversion :

( )

L’équation bilan par rapport à l’acide (cas d’un réacteur à cuve agité fermé) est exprimée par :

( ) 

 ( ) Equation bilan de matière

( ) ( )

Bilan d’énergie par rapport à chaque kilogramme de l’huile :

Cas du réacteur adiabatique :

Bilan entre t0=0 et t :

∆ ( )
∆ ( ) ( )  ( )
Nombre initial de moles d’acide contenu dans chaque kg d’huile : N0= mA/MA=0.156/60
Ainsi :
Equation bilan thermique cas du réacteur adiabatique

Cas du réacteur non-adiabatique :

∆ ( ) ( ) ( )

Ainsi :
Equation bilan thermique cas du réacteur non-adiabatique

28
 Ecole Mohammadia d’Ingénieurs- GPI 2

Résolution Numérique: Méthode de Runge-Kutta :

Réacteur Adiabatique :
Conditions initiales : X0=0 ; T0=613, t0=0 et incrément ∆X=0.02
Pour chaque itération n+1 :

∆ ( ( ))
∆ ∆
∆ ( ( ))
∆ ∆
∆ ( ( ))
∆ ( ∆ ( ∆ ))
( ) ( )
∆

Créations des programmes sur Matlab :

Fonction g1

function g1=g1(X,T)
R=1.987; % en cal/mol/K
g1=60./((1-X).*exp(35-44500./(R*T)));
end

Fonction de température T(X)

function T1=T1(X)
T1=613-65*X;
end

Programme de résolution :

clear all
clc

size=35; % Nombre d'itérations

t=zeros(1,size); % rapport entre volume du réacteur et débit d'alimentation
en B: V/FBf
X=zeros(1,size);
T=zeros(1,size);
T(1)=613;
DX=0.02; % Incrément pour V/FBf

for i=2:size
k0=DX*g1(X(i-1),T1(X(i-1)));
k1=DX*g1(X(i-1)+DX/2,T1(X(i-1)+DX/2));
k2=DX*g1(X(i-1)+DX/2,T1(X(i-1)+DX/2));
k3=DX*g1(X(i-1)+DX,T1(X(i-1)+DX));

t(i)=t(i-1)+(1/6)*(k0+(2*k1)+(2*k2)+k3);
X(i)=X(i-1)+DX;
T(i)=613-65*X(i);

29
 Ecole Mohammadia d’Ingénieurs- GPI 2

end

plot(t/60,X,'-x')
xlabel('Temps en (min)')
ylabel('Taux de Conversion')
grid
figure
plot(t/60,T,'-*')
xlabel('Temps en (min)')
ylabel('Température (K)')
grid

Résultats obtenus :

Réacteur non-Adiabatique :
Conditions initiales : X0=0 ; T0=613, t0=0 et incrément ∆X=0.02
Pour chaque itération n+1 :

∆ ( ( ))
∆ ∆
∆ ( ( ))
∆ ∆
∆ ( ( ))
∆ ( ∆ ( ∆ ))
( ) ( )
∆

30
 Ecole Mohammadia d’Ingénieurs- GPI 2

Création d’une nouvelle fonction de température :

function T2=T2(X,t)
T2=613-65*X+0.0922*t;
end

Programme de résolution :
size=20; % Nombre d'itérations
t=zeros(1,size); % rapport entre volume du réacteur et débit d'alimentation
en B: V/FBf
X=zeros(1,size);
T=zeros(1,size);
T(1)=613;
DX=0.02; % Incrément pour V/FBf

for i=2:size
k0=DX*g1(X(i-1),T2(X(i-1),t(i-1)));
k1=DX*g1(X(i-1)+DX/2,T2(X(i-1)+DX/2,t(i-1)+k0/2));
k2=DX*g1(X(i-1)+DX/2,T2(X(i-1)+DX/2,t(i-1)+k1/2));
k3=DX*g1(X(i-1)+DX,T2(X(i-1)+DX,t(i-1)+k2));

t(i)=t(i-1)+(1/6)*(k0+(2*k1)+(2*k2)+k3);
X(i)=X(i-1)+DX;
T(i)=613-65*X(i)+0.0922*t(i);
end

Résultats obtenus :

31
 Ecole Mohammadia d’Ingénieurs- GPI 2

 Modélisation d’un collecteur Solaire Cylindro-Parabolique (Solar PTC )

Figure 1: Simplified schematic of PTC collector

The parabolic trough collector is a one-dimensional parabola that focuses the parallel sun rays into a focal
line (Fig. 1.). This collector is equipped with a receiver tube positioned in the focal line, and it absorbs the
concentrated solar rays. A significant reduction of heat dissipated from the receiver tube is obtained by
covering the receiver by an evacuated glass tube called glass cover. Heat transmitted through the receiver
tubes walls is absorbed by a heating medium (synthetic oil, water, molten salts…etc.). The reflector
consists of a highly reflecting parabolic mirror. The optical standards specifications of the reflector and
receiver are provided in Table below.
Table1: optical standards specifications of the reflector and receiver

The parabolic through collector is used in order to collect and convert the solar energy into thermal
energy, the useful energy produced by the system could be calculated by:
( ( )) ( )

32
 Ecole Mohammadia d’Ingénieurs- GPI 2

Where Aap, Ar, FR and UL are respectively the aperture area, the receiver area and the solar collector heat
loss coefficient between ambient and receiver. The receiver efficiency and the heat absorbed by the
receiver are given by equations 8 and 9:
( )

( )

The aperture area is calculated in function of the geometrical characteristics of the parabolic trough
collector:
( ) ( )

The heat removal factor is defined as:

( ( )) ( )
(

The solar collector heat loss coefficient between ambient and receiver is expressed with the equation:
( )
( )

Where hr,ca is the radiation heat transfer coefficient between ambient and the cover expressed with:
( )( ) ( )

And hr,cr is the radiation heat transfer coefficient between the receiver and the cover, it is defined by:
( )( )
( )
( )

Finally hc,ca represents the convective heat transfer coefficient between the cover and the air of ambient, it
is calculated in function of Nusselt number using the equation below:
( ) ( )

The overall heat coefficient from the surroundings to the heating fluid flowing inside the receiver is
defined as following:
( ( )) ( )

Where hc,r,in is the convective heat transfer coefficient between the heating medium and the internal
receiver wall, it is defined by:
( ) ( )

The temperature of the cover is calculated using the equation:

( )
( )
( )

The thermal energy which is transferred to the heat transfer fluid is also equal to:
( )

33
 Ecole Mohammadia d’Ingénieurs- GPI 2

Nusselt number for air could be calculated by:

For 0.1 < Re < 1000:

For 1000 < Re < 50 000:

And Nusselt number for heating medium flowing inside the receiver tube could be evaluated by:

The geometrical specifications of each single PTC collector are given the following table:
Table2: Standard Geometrical specifications of PTC

 Work to do using:
 Modeling based on physical and thermal phenomena
 Statistical Modeling

The PTC collector uses synthetic oil as heating medium, the oil enters the collector at 150°C with a
constant volume flow rate equals to 1.2.10-3m3/s (1.04kg/s). Wind velocity is 4m/s and ambient air
temperature is 20°C.
1) We assume that the collector operates under steady state conditions. Calculate the overall collector
efficiency for a constant solar irradiation (DNI): DNI=900 W/m2;
2) Represents graphically the change of overall thermal efficiency in relation with the solar
irradiation (within a range of 500-1300 W/m2);
3) Represents in the same figure (3d) the change of overall thermal efficiency versus the receiver
diameter and the collector length. (Length: 4-15 m , inner receiver diameter: 0.015-0.04 m).
4) Determine the optimal geometrical specifications of the collector within ranges: Length: 10-16m,
Inner receiver diameter: 0.015-0.04m, Collector width: 2.5-5m. The heating medium temperature at
the outlet must not exceed 200°C. The 24hours solar irradiation is given by: DNI=[0,0,0,0,0,
0,400,700,800,900,1000,1200,1200,1000,900,800,700,500,400,0,0,0,0,0].

Additional data: Specifications of synthetic oil used as heating medium

 Cp=1589+3.206*T-5.462*(10^(-4))*(T.^2); % Specific heat capacity(J/kg/K)
 CinVisc=0.0045*exp(-0.0158*T); % Viscosité cinématic (m2/s)
 desity=980-0.759*T+3.214*(10^(-6))*(T.^2);% desity (kg/m3)
 viscoo=desity.*CinVisc; % Dynamic viscosity of oil (kg/(m.S))
 ko=0.125-1.391*(10^(-4))*T-4.686*(10^(-8))*(T.^2); % Conductivity of oil W/K/m
Reference:
John [Link] and william [Link], Solar Engineering of Thermal Processes, Second edi. John Wiley & Sons.
J. A. Duffie and W. A. Beckman, Solar Engineering of Thermal Processes. 2006.
S. A. Kalogirou, Solar Energy Engineering Processes and Systems, Second Ed. Elsevier, 2013.
T. L. Bergman, A. S. Lavigne, F. P. Incropera, and D. P. Dewitt, Fundamentals of Heat and Mass Transfer, Seventh ed. 2011.

34
 Ecole Mohammadia d’Ingénieurs- GPI 2

Experimental data (Data to be used for statistical modeling):

(See: Excel File: PTC DATA)

 Resolution algorithm (to be discussed with students)

 Code using Matlab:

Calculation script (function):

function
PTCfct=PTCfct(AbsReceiver,DNI,W,dce,EmitCover,EmitReceiver,ReflMirror,dre,Tra
nsGlass,Te,Tamb,M)

L=12; % length (m)

dri=dre-0.01; % Inner Receiver diameter (m)
kmetal=20; % W./m./K
IntrFact=0.93; % Intercept factor
IncAngl=1; % Incidence angle modifier
%--------------------------------------------------------------------------
% Données opératoires et thermiques/thermo...
vair=10; % Ambient air Velocity (m/s)
roair=1.176; % air density (kg/m3)
kair=0.0249; % en (W/m/K)
Viscair=1.823.*10^(-5); % kg/m/s
Reair=roair.*vair.*dce./Viscair; % Nombre de reynolds pour l'air
% Optical effeciency calculation
EffOptical=EmitCover*ReflMirror*IntrFact*TransGlass*AbsReceiver*IncAngl;
S=DNI.*EffOptical;
% Estimation des températures
Tr0=200; % En (°C)
Ts0=180; % En (°C)
Tc0=200; % En (°C)
%--------------------------------------------------------------------------

Error1=1000; % Error for T Receiver calculation

while Error1>1 % While loop for Receiver temperature
Tr=Tr0;
Error2=1000; %Error for T out calculation
while Error2>1 % While loop for outlet temperature
Ts=Ts0;
Error3=1000; %Error for T cover calculation
while Error3>1 % While loop for cover temperature
Tc=Tc0;
Tmoy=((Ts+Te)/2)+273;
% Thermal and physical specifications of heating medium
Cp=1589+3.206*Tmoy-5.462*(10^(-4))*(Tmoy.^2); % Specific heat
capacity(J/kg/K)
CinVisc=0.0045*exp(-0.0158*Tmoy); % Viscosité cinématic (m2/s)
desity=980-0.759*Tmoy+3.214*(10^(-6))*(Tmoy.^2);% desity
(kg/m3)
viscoo=desity.*CinVisc; % Dynamic viscosity of oil (kg/(m.S))

35
 Ecole Mohammadia d’Ingénieurs- GPI 2

ko=0.125-1.391*(10^(-4))*Tmoy-4.686*(10^(-8))*(Tmoy.^2); %
Conductivity of oil W/K/m
vf=4.*M./(desity.*pi.*dri.^2);
% Calcul de hr,ca
hrca=EmitCover.*5.67*(10^(-
8)).*(2*273+Tc+Tamb).*((Tc+273).^2+(Tamb+273).^2); % Radiation between
ambient and cover (W/m2/K)
hrcr=5.67*10^(-
8).*(2*273+Tc+Tr).*((Tc+273).^2+(Tr+273).^2)./((1./EmitReceiver)+(dre./dce).*
((1./EmitCover)-1)); % Radiation between Receiver and cover

XX=Reair;

if XX<1000
Nuair=0.4+0.54*Reair.^0.52;
else
Nuair=0.3*Reair.^0.6;
end
hcca=Nuair.*kair./dce; % Convection heat transfer coefficient
(W/m2/K)
Ref=desity.*vf.*dri./viscoo;
Prf=Cp.*viscoo./ko;
Nuf=0.023*(Ref.^0.8).*(Prf.^0.4);
hfi=Nuf.*ko./dri; % Convection fluide circulant à l'intérieur
du tube absorbeur (W/m2/K)
UL=((dre./dce).*(1./(hcca+hrca))+(1./hrcr))^(-1); % Exprimée en
W/m2/K (Coefficient Receiver-Ambiant)
U0=(
(1./UL)+((dre./dri).*(1./hfi))+((dre./(2*kmetal)).*log(dre./dri))).^(-1); %
Exprimée en W/m2/K (Coefficient Fluide- Ambiant)
Ar=pi.*dre.*L;
FR=(M.*Cp./(Ar.*UL)).*(1-exp(-Ar.*U0./(M.*Cp))); % Removal
factor
Ac=pi.*dce.*L;

Tc0=(hrcr.*(Tr+273)+(Ac./Ar).*(hcca+hrca).*(Tamb+273))./(hrcr+(Ac./Ar).*(hcca
+hrca))-273;
Error3=abs(Tc0-Tc);

end
Aa=W.*L;
Qu=Aa.*FR.*(S-(Ar./Aa).*UL.*(Tr-Tamb));
Ts0=Te+Qu./(M.*Cp); % Nouveau calcul de T sortie
Error2=abs(Ts0-Ts);
Ts=Ts0;
end
Flux=Qu./Ar;
Rmetal=((dre./(2*kmetal)).*log(dre./dri));
Tr0=Tmoy-273+Flux.*((1./hfi)+Rmetal);
Error1=abs(Tr0-Tr);

end

Eff=Qu/(Aa*DNI);

36
 Ecole Mohammadia d’Ingénieurs- GPI 2

PTCfct=[Ts,Eff];
end

Results Plotting Script

% Variables
AbsReceiver=xlsread('PTCData','B:B')/100;
DNI=xlsread('PTCData','C:C');
W=xlsread('PTCData','D:D');
dce=xlsread('PTCData','E:E');
EmitCover=xlsread('PTCData','F:F')/100;
EmitReceiver=xlsread('PTCData','G:G')/100;
ReflMirror=xlsread('PTCData','H:H')/100;
dre=xlsread('PTCData','I:I');
TransGlass=xlsread('PTCData','J:J')/100;
Te=xlsread('PTCData','K:K');
Tamb=xlsread('PTCData','L:L');
M=xlsread('PTCData','M:M');

% les données mésurées-Output

Ts_exp=xlsread('PTCData','O:O');
Effeciency_exp=xlsread('PTCData','P:P');

% Utilisation du modèle de calcul -----------------------------------------

Size=size(M);
N=Size(1);
Ts=zeros(N,1);
Effeciency=zeros(N,1);
Cpp=zeros(N,1);
for i=1:N

Results=PTCfct(AbsReceiver(i),DNI(i),W(i),dce(i),EmitCover(i),EmitReceiver(i)
,ReflMirror(i),dre(i),TransGlass(i),Te(i),Tamb(i),M(i));
Ts(i)=Results(1);
Effeciency(i)=Results(2);
end

% Représentation des résultats

figure
plot(Ts_exp,Ts,'*',70:1:150,70:1:150,'-r')
xlabel('Ts Exper. (°C)')
ylabel('Ts Modèle')
grid

figure
plot(100*Effeciency_exp,100*Effeciency,'*',0:1:70,0:1:70,'-r')
xlabel('Effeciency Exper. (%)')
ylabel('Effeciency Modèle (%)')
grid

37
 Ecole Mohammadia d’Ingénieurs- GPI 2

Experimental validation

Figure 2: Experimental validation- Outlet temperature

Figure 3: Experimental validation-Overall efficiency

38
 Ecole Mohammadia d’Ingénieurs- GPI 2

6-1- Dimensionnement d’un évaporateur

Un évaporateur à tubes horizontaux et à film tombant, utilisant l’eau comprimée comme fluide
caloporteur pour évaporer, à basse pression, une portion de l’eau douce d’alimentation. Les propriétés
thermiques des deux fluides sont résumées dans le tableau suivant :

Propriétés Eau chaude sous-Pression Eau d’alimentation

Température d’Entrée (°C) 100 85
Température de Sortie (°C) 88 85
Chaleur latente d’évaporation 2295 kJ/kg

Les tubes de l’évaporateur ont une longueur de 1 m, le débit de l’eau chaude comprimée est de 100 kg/s
et le taux d’évaporation de l’eau d’alimentation est de 40%. L’évaporateur est estimé adiabatique, et la
résistance d’encrassement totale est de l’ordre de 0.4 (m2K)/kW.

1- Calculer la surface d’échange de l’évaporateur et le coefficient d’échange thermique global pour les
diamètres des tubes suivants :

Diamètre intérieur (mm) Epaisseur (mm) Dénomination (pouce) Matériau

26.9 2.6 ¾ ’’ Acier-10% Ni
42.4 3.2 1’’ 1/4 Acier-10% Ni
48.3 3.2 1’’ 1/2 Acier-10% Ni

2- Représenter graphiquement l’effet de la température de l’eau chaude comprimée et le diamètre des

tubes sur les coefficients d’échange thermique et la surface d’échange thermique.

Marge de variation de température : 95-120°C. Marge de variation du diamètre : 10-40mm.

3- Trouver le diamètre et la température d’évaporation optimisés en utilisant la méthode de l’algorithme

génétique.

39
 Ecole Mohammadia d’Ingénieurs- GPI 2

Données :

Propriétés thermiques de l’eau liquide:

Cp_l=4135.1+2.5*T-3*(10^(-2))*(T.^2)+(10^(-4))*(T.^3);
CinVisc_l =(10^(-5))*(T.^(-0.838));
desity_l =1002.5-0.1714*T-0.0026*(T.^2);
viscoo_l =desity_l.*CinVisc_l;
ko_l =0.5685+0.0016*T-6*(10^(-6))*(T.^2);

Corrélations pour les Coefficient thermique:

Dans le cas de l’eau chaude liquide circulant à l’intérieur des tubes lisses ; le coefficient d’échange
thermique par convection est estimé en utilisant la corrélation développée par Kern [1]:
 Pour des valeurs de Re<2000 :
hi⁄ dG µ µ
G ( i ⁄µ) ⁄
( ⁄ ) ⁄
(L⁄d ) ⁄
( ⁄µ )
i

 Pour des valeurs de 2000⫹Re<10000:

⁄3 −
dG (L⁄d ) 3
hi⁄ ( i ⁄µ)
i
G ( (di G⁄ )( )
−
µ) cµ µ
( ⁄k)3 ( ⁄µw )

 Pour les valeurs de Reynolds supérieures à 10000:

hi⁄ µ ⁄ dG µ
G ( ⁄ ) ( i ⁄µ) ( ⁄µ )

Dans le cas d’ébullition d’un film tombant liquide circulant à travers la surface extérieure d’un faisceau
tubulaire horizontale, le coefficient d’échange thermique par convection peut être estimé en utilisant la
corrélation de Mu et Shen [2]:
µ
h (υ ⁄ ) ( Γ⁄µL ) ( L ⁄ ) 7
(B.15)

Référence:

[1] D.Q. Kern, Process Heat Transfer, International Student Edition, 1950.
[2] Xingsen Mu , Shengqiang Shen , Yong Yang & Xiaohua Liu (2012) Experimental study of falling film evaporation
heat transfer coefficient on horizontal tube, Desalination and Water Treatment, 50:1-3, 310-316, DOI:
10.1080/19443994.2012.719734.

40
 Ecole Mohammadia d’Ingénieurs- GPI 2

Algorithme de calcul :

Input Known Parameters

Calculate LMTD

Estimate U=U0
Replace U estimated by U
calculated

Calculate H.T. Area

Calculate Number of Tubes

NO
U =estimated U ?

Calculate the Convective

YES
Heat Transfer Coefficients

End
Calculate the new value of U

r
r ( ⁄rin )
U R
hin rin h

Code de Dimensionnement :

% Etape 1: Données du problème:

Teh=100+273; % Température d'entrée de l'eau chaude en K

Tsh=88+273; % Température de sortie de l'eau chaude en K
Teau=85+273; % Température d'évaporation de l'eau en K
Hev=2400; % Chaleur latente d'évaporation de l'eau (kJ/kg)
Mh=100; % Débit de l'eau thermique (kg/s)
Xev=0.4; % Fraction de l'eau évaporée
Rf=0.4; %%% fouling resistance in (m2.°C/kW)

Kmetal=26/1000; % Conductivité thermique du métal en W/m/K

di=0.0269; % Diamètre interne des tubes en (m)

41
 Ecole Mohammadia d’Ingénieurs- GPI 2

e=0.0026; % Epaisseur du tube en (m)

do=di+2*e; % Diamètre externe en (m)
L=1; % Longueur du tube en (m)
section=pi.*(di.^2)./4; % Section de chaque tube en m2

% Etape 2: Propriétés thermiques,Bilans thermiques et calcul de LMTD

Tmh=(Teh+Tsh)/2;
Cphmoy=(4135.1+2.5*Tmh-3*(10^(-2))*(Tmh.^2)+(10^(-4))*(Tmh.^3))/1000; %
Capacité thermique moyenne de l'eau chaude exprimée en kJ/kg/K
CinVisch=(10^(-5))*(Tmh.^(-0.838)); % Viscosité cinématique moyenne de l'eau
chaude exprimée en m2/s
Densityh=1002.5-0.1714*Tmh-0.0026*(Tmh.^2); % densité moyenne de l'eau chaude
exprimée en kg/m3
viscooh=Densityh.*CinVisch; % Viscosity dynamique moyenne de l'eau chaude
exprimée en kg/m/s
Kh=(0.5685+0.0016*Tmh-6*(10^(-6))*(Tmh.^2))/1000; % Conductivité moyenne de
l'eau chaude en kW/m/K

Cpw=(4135.1+2.5*Teau-3*(10^(-2))*(Teau.^2)+(10^(-4))*(Teau.^3))/1000; %
Capacité thermique moyenne de l'eau exprimée en kJ/kg/K
CinViscw=(10^(-5))*(Teau.^(-0.838)); % Viscosité cinématique moyenne de l'eau
exprimée en m2/s
Densityw=1002.5-0.1714*Teau-0.0026*(Teau.^2); % densité moyenne de l'eau
exprimée en kg/m3
viscoow=Densityw.*CinViscw; % Viscosity dynamique moyenne de l'eau exprimée
en kg/m/s
Kw=(0.5685+0.0016*Teau-6*(10^(-6))*(Teau.^2))/1000; % Conductivité moyenne de
l'eau en kW/m/K

Q=Mh.*Cphmoy.*(Teh-Tsh)/1000; % Puissance échangée en kJ/s ou kW

Mvap=Q./Hev; % Débit de la vapeur produite (kg/s)
Meau=Mvap./Xev; % Débit de l'eau d'alimentation en kg/s
LMTD=((Teh-Teau)-(Tsh-Teau))./log((Teh-Teau)./(Tsh-Teau));

% Etape 3: Boucle de calcul du coefficient d'échange thermique

U0=0.9; % Première estimation du coefficient d'échange thermique
Error=1000; % Initialisation de l'erreur

while Error>0.0001
U1=U0;
A=Q./(U1.*LMTD); % Surface d'échange thermique (m2)
Ntubes=A./(do.*L.*pi); % Nombre des tubes
Mht=Mh./Ntubes; % Débit massique de l'huile par tube

% convection à l'extérieur des tubes

gama=((Meau./Ntubes)./L); % Débit massique par longueur de tube kg/m
Y=0.0532*((4*gama./viscoow).^0.21).*(((Cpw.*viscoow)./Kw).^0.731).*(exp(-
0.02283*Xev));
ho=Y.*(((CinViscw.^2)./(9.18*(Kw.^3))).^(-1/3));
Rboiling=1./ho; % Résistance thermique (m2 K/W )

42
 Ecole Mohammadia d’Ingénieurs- GPI 2

% Convection à l'intérieur des tubes

Re=(4*Mh)./((3.14*di).*Ntubes.*viscoow); % nombre de reynolds

G=(4*Mh)./(3.14*di.^2); % Flux massique kg/s/m2

if Re(1)<2000
hi=(1.86*Cphmoy.*G.*((di.*G./viscooh).^(-
2/3)).*((Cphmoy.*viscooh./Kh).^(-2/3)).*((L./di).^(-
1/3)).*((viscooh./viscoow).^0.14))/1000; %kW
elseif Re(1)>2000
if Re(1) <10000
hi=(0.116*Cphmoy.*G.*((((di.*G./viscooh).^(2/3))-
125)./(di.*G./viscooh)).*((1+((L./di).^(-
2/3)))./(((Cphmoy.*viscooh./Kh).^(2/3)).*((viscooh./viscoow).^0.14))))/1000;
%kW
else
hi=(0.023*Cphmoy.*G.*((di.*G./viscooh).^(-
0.2)).*((Cphmoy.*viscooh./Kh).^(-2/3)).*(1+((L./di).^(-
0.7))).*((viscooh./viscoow).^0.14))/1000; % kW
end
end

Rcooling=(do./di)./hi; % Résistance thermique (m2 K/W )

Rmetal=(do./2).*log(do/di)./Kmetal; % Résistance du métal

U0=1./(Rcooling+Rboiling+Rf+Rmetal);

Error=abs(U0-U1);
end

Code de Dimensionnement et analyse thermique : Effet de Température et diamètre

% Etape 1: Données du problème:

Teh=linspace(95+273,130+273,20); % Température d'entrée de l'huile thermique

en K
di=linspace(0.01,0.04,20); % Diamètre interne des tubes en (m)
[Teh,di]=meshgrid(Teh,di);
Tsh=88+273; % Température de sortie de l'huile thermique en K
Teau=85+273; % Température d'évaporation de l'eau en K
Hev=2400; % Chaleur latente d'évaporation de l'eau (kJ/kg)
Mh=100*ones(20,20); % Débit de l'eau thermique (kg/s)
Xev=0.4; % Fraction de l'eau évaporée
Rf=0.4; %%% fouling resistance in (m2.°C/kW)

Kmetal=26/1000; % Conductivité thermique du métal en W/m/K

e=0.0026; % Epaisseur du tube en (m)
do=di+2*e; % Diamètre externe en (m)
L=1; % Longueur du tube en (m)
section=pi.*(di.^2)./4; % Section de chaque tube en m2

43
 Ecole Mohammadia d’Ingénieurs- GPI 2

% Etape 2: Propriétés thermiques,Bilans thermiques et calcul de LMTD

Tmh=(Teh+Tsh)/2;
Cphmoy=(4135.1+2.5*Tmh-3*(10^(-2))*(Tmh.^2)+(10^(-4))*(Tmh.^3))/1000; %
Capacité thermique moyenne de l'eau chaude exprimée en kJ/kg/K
CinVisch=(10^(-5))*(Tmh.^(-0.838)); % Viscosité cinématique moyenne de l'eau
chaude exprimée en m2/s
Densityh=1002.5-0.1714*Tmh-0.0026*(Tmh.^2); % densité moyenne de l'eau chaude
exprimée en kg/m3
viscooh=Densityh.*CinVisch; % Viscosity dynamique moyenne de l'eau chaude
exprimée en kg/m/s
Kh=(0.5685+0.0016*Tmh-6*(10^(-6))*(Tmh.^2))/1000; % Conductivité moyenne de
l'eau chaude en kW/m/K

Cpw=(4135.1+2.5*Teau-3*(10^(-2))*(Teau.^2)+(10^(-4))*(Teau.^3))/1000; %
Capacité thermique moyenne de l'eau exprimée en kJ/kg/K
CinViscw=(10^(-5))*(Teau.^(-0.838)); % Viscosité cinématique moyenne de l'eau
exprimée en m2/s
Densityw=1002.5-0.1714*Teau-0.0026*(Teau.^2); % densité moyenne de l'eau
exprimée en kg/m3
viscoow=Densityw.*CinViscw; % Viscosity dynamique moyenne de l'eau exprimée
en kg/m/s
Kw=(0.5685+0.0016*Teau-6*(10^(-6))*(Teau.^2))/1000; % Conductivité moyenne de
l'eau en kW/m/K

Q=Mh.*Cphmoy.*(Teh-Tsh)/1000; % Puissance échangée en kJ/s ou kW

Mvap=Q./Hev; % Débit de la vapeur produite (kg/s)
Meau=Mvap./Xev; % Débit de l'eau d'alimentation en kg/s
LMTD=((Teh-Teau)-(Tsh-Teau))./log((Teh-Teau)./(Tsh-Teau));

% Etape 3: Boucle de calcul du coefficient d'échange thermique

U0=0.9*ones(20,20); % Première estimation du coefficient d'échange thermique
Error=1000*ones(20,20); % Initialisation de l'erreur

U1=zeros(20,20);
A=zeros(20,20);
Ntubes=zeros(20,20);
Mht=zeros(20,20);
gama=zeros(20,20);
ho=zeros(20,20);
Rboiling=zeros(20,20);
Re=zeros(20,20);
Y=zeros(20,20);
hi=zeros(20,20);
G=zeros(20,20);
Rcooling=zeros(20,20);
Rmetal=zeros(20,20);

for i=1:20
for j=1:20

while Error(i,j)>0.0001
U1(i,j)=U0(i,j);

44
 Ecole Mohammadia d’Ingénieurs- GPI 2

A(i,j)=Q(i,j)./(U1(i,j).*LMTD(i,j)); % Surface d'échange thermique (m2)

Ntubes(i,j)=A(i,j)./(do(i,j).*L.*pi); % Nombre des tubes
Mht(i,j)=Mh(i,j)./Ntubes(i,j); % Débit massique de l'huile par tube
% convection à l'extérieur des tubes
gama(i,j)=((Meau(i,j)./Ntubes(i,j))./L); % Débit massique par longueur de
tube kg/m

Y(i,j)=0.0532*((4*gama(i,j)./viscoow).^0.21).*(((Cpw.*viscoow)./Kw).^0.731).*
(exp(-0.02283*Xev));
ho(i,j)=Y(i,j).*(((CinViscw.^2)./(9.18*(Kw.^3))).^(-1/3));
Rboiling(i,j)=1./ho(i,j); % Résistance thermique (m2 K/W )

% Convection à l'intérieur des tubes

Re(i,j)=(4*Mh(i,j))./((3.14*di(i,j)).*Ntubes(i,j).*viscooh(i,j)); % nombre
de reynolds
G(i,j)=(4*Mh(i,j))./(3.14*di(i,j).^2); % Flux massique kg/s/m2

if Re(i,j) <2000
hi(i,j)=(1.86*Cphmoy(i,j).*G(i,j).*((di(i,j).*G(i,j)./viscooh(i,j)).^(-
2/3)).*((Cphmoy(i,j).*viscooh(i,j)./Kh(i,j)).^(-2/3)).*((L./di(i,j)).^(-
1/3)).*((viscooh(i,j)./viscoow(i,j)).^0.14))/1000; %kW
elseif Re(i,j) >2000
if Re(i,j)<10000

hi(i,j)=(0.116*Cphmoy(i,j).*G(i,j).*((((di(i,j).*G(i,j)./viscooh(i,j)).^(2/3)
)-125)./(di(i,j).*G(i,j)./viscooh(i,j))).*((1+((L./di(i,j)).^(-
2/3)))./(((Cphmoy(i,j).*viscooh(i,j)./Kh(i,j)).^(2/3)).*((viscooh(i,j)./visco
ow).^0.14))))/1000; %kW
else

hi(i,j)=(0.023*Cphmoy(i,j).*G(i,j).*((di(i,j).*G(i,j)./viscooh(i,j)).^(-
0.2)).*((Cphmoy(i,j).*viscooh(i,j)./Kh(i,j)).^(-2/3)).*(1+((L./di(i,j)).^(-
0.7))).*((viscooh(i,j)./viscoow).^0.14))/1000; % kW
end
end
Rcooling(i,j)=(do(i,j)./di(i,j))./hi(i,j); % Résistance thermique (m2 K/W )
Rmetal(i,j)=(do(i,j)./2).*log(do(i,j)/di(i,j))./Kmetal; % Résistance du métal
U0(i,j)=1./(Rcooling(i,j)+Rboiling(i,j)+Rf+Rmetal(i,j));
Error(i,j)=abs(U0(i,j)-U1(i,j));
end
end
end
figure
mesh(Teh-273,1000*di,U1)
xlabel('Température The (°C)')
ylabel('Diamètre interne du tube(mm)')
title('Coefficient d echange global(kW/m2/K)')
figure
mesh(Teh-273,1000*di,hi)
xlabel('Température The (°C)')
ylabel('Diamètre interne du tube(mm)')
title('Coefficient d echange hi (kW/m2/K)')
figure
mesh(Teh-273,1000*di,ho)
xlabel('Température The (°C)')

45
 Ecole Mohammadia d’Ingénieurs- GPI 2

ylabel('Diamètre interne du tube(mm)')

title('Coefficient d echange ho (kW/m2/K)')

Résultats:

46