Discussion :

[Mani8]

Bonjour,

Je cherche à implémenter dans Matlab un algorithme permettant de "résoudre" un réseau de tuyauterie (c'est à dire déterminer le débit Q et la direction du flux dans chaque tuyau ainsi que la hauteur piézométrique H à chaque intersection).
Le résultat obtenu est bien celui que j'attends mais pour une raison que je ne parviens pas à identifier la vitesse de convergence de l'algorithme est très lente (en témoignent les graphes de la déviation relative maximale et du résidu maximal).

Quelqu'un aurait-il un conseil pour accélérer la convergence ? Mon problème vient-il d'une fonction que j'aurais mal utilisé ?

Le code:

Code matlab :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
clear
clear all
close all
 
% Fluid and gravity
rho = 1e3;
mu  = 1e-3;
g   = 9.81;
 
% Colebrook formula
colebrook = @(Re,rr,la) 1/sqrt(la)+2*log10(rr/3.71+2.51/Re/sqrt(la));
 
% Deposit
h1 = 200;
h2 = 200;
 
% Pipe data
L1 = 600;
L2 = 400;
L3 = 900;
L4 = 400;
L6 = 350;
L7 = 600;
L9 = 300;
L10 = 500;
L11 = 300;
L12 = 350;
 
D1 = 0.5;
D2 = 0.5;
D3 = 0.5;
D4 = 0.5;
D6 = 0.5;
D7 = 0.5;
D9 = 0.3;
D10 = 0.3;
D11 = 0.3;
D12 = 0.3;
 
v_L  = [L1 L2 L3 L4 L6 L7 L9 L10 L11 L12];
v_D  = [D1 D2 D3 D4 D6 D7 D9 D10 D11 D12];
 
v_r = 1e-3*ones(1,size(v_L,2));
v_rr = v_r./v_D;
 
% Initial lambda values. Starting with Reynolds infinity.
v_la = 0.25./(log10(v_rr/3.71)).^2;
 
% Initial conditions
v_Q = 0.1*ones(1,size(v_L,2));
v_H = [10 h1 h2 10 10 10 10];
 
%Connectivity matrix
    % 1  2  3  4  6  7  8 
E = [-1  1  0  0  0  0  0; %1
     -1  0  1  0  0  0  0; %2
      1  0  0 -1  0  0  0; %3
      0  1  0 -1  0  0  0; %4
      0  0  1  0 -1  0  0; %6
      0  0  0  1 -1  0  0; %7
      0  0  1  0  0 -1  0; %9
      0  0  0  0  1 -1  0; %10
      0  0  0  0  1  0 -1; %11
      0  0  0  0  0  1 -1];%12
 
 %Declare known C values 
 x=0;
 v_C=[-0.2 x x -0.3 0 0 -0.5];
 y = [0 1 1 0 0 0 0]; %logic matrix: 0 when C is known, 1 when H is known
 
 %%....Main Loop....%%
 
 %Iteration limits
 maxiter = 1000;
 tol = 1e-3;
 inc = 1;
 res = 1;
 iter = 0;
 
 %Creation of solution storage matrices
 Q_results = zeros(size(v_Q,2),maxiter);
 H_results = zeros(size(v_H,2),maxiter);
 inc_results = zeros(1,maxiter);
 res_results = zeros(1,maxiter);
 v_sol = [v_Q v_H(1,1) v_H(1,4:size(v_H,2))];
 
 %...........MAIN LOOP.............%
 while (inc>tol) && (res>tol) && (iter<maxiter)
 
 iter = iter+1;
 
 %Storage of previous iteration results
 Q_results(:,iter)=v_Q';
 H_results(:,iter)=v_H';
 v_sol_old = v_sol;
 
 %Calculate Reynolds numbers from Q
 v_vel = 4*v_Q/pi./v_D.^2;
 v_Re  = abs(rho*v_vel.*v_D/mu);
 
 %Get the new lambdas
 for i = 1:size(v_Q,2)
    v_la(i) = fzero(@(x)colebrook(v_Re(i),v_rr(i),x),v_la(i));
 end
 
 %Obtain matrix M
 v_F = 8/pi^2/g./v_D.^4.*v_la.*v_L./v_D;
 v_K = v_F.*abs(v_Q);
 K = diag(v_K);
 M = E'/K*E;
 
 %Obtain A and b
 A=M*diag(1-y)-diag(y);
 b=-M*(v_H.*y)'+(v_C.*(1-y))';
 
 %Solve the system
 X=A\b;
 
 %Extract the new H and C values 
 v_H=X';
 v_H(2)=h1;
 v_H(3)=h2;
 v_C(2)=X(2);
 v_C(3)=X(3);
 
 %Re-evaluate Q
 v_Q=(K\E*v_H')';
 
 %Update solution matrix
 v_sol = [v_Q v_H(1,1) v_H(1,4:size(v_H,2))];
 
 %Relative increment of solution
 inc = norm((v_sol-v_sol_old)./v_sol,inf);
 %inc = norm((v_Q-v_sol_old(1:10))./v_Q,inf);
 inc_results(iter) = inc;
 
 %Residual of system
 res = norm(v_sol-v_sol_old,inf);
 %res = norm(v_Q-v_sol_old(1:10),inf);
 res_results(iter) = res;
 
 end
 
 %Plots
 inc_results(iter+1:maxiter) = [];
 res_results(iter+1:maxiter) = [];
 
 figure(1)
 subplot(211)
 semilogy(1:iter,inc_results)
 box on,grid on
 xlabel('iteration')
 ylabel('relative increment')
 
 subplot(212)
 semilogy(1:iter,res_results)
 box on,grid on
 xlabel('iteration')
 ylabel('residue')

Le système à résoudre:


Détail de l'algorithme:
1- Première estimation des valeurs Q et H + estimation arbitraire de l'orientation des flux (correspond aux flèches bleues sur le schéma). Ces orientations sont consignées dans une matrice de connectivité c.a.d une matrice ayant autant de colonnes qu'il y a d'embranchements et autant de lignes qu'il y a de tuyaux. Les éléments de cette matrice valent 1 là d'où le flux part pour un tuyau donné, -1 là où il arrive et 0 partout ailleurs. Cette matrice s’appellera E.

2- Calcul de la vitesse d'écoulement v dans chaque tuyau à partir de Q=(pi*D^2/4)*v

3- Calcul du nombre de Reynolds pour chaque tuyau

4- Utilisation de la formule de Colebrook pour déterminer le coefficient de Darcy (pour la première utilisation de la formule, j'utilise une valeur du coefficient de Darcy correspondant à un nombre de Reynolds tendant vers l'infini c.a.d.
lambda=0.25/( log(eps/3.71)^2 ) )


5- Pour chaque tuyau, je calcule la valeur K=lambda/(2*g*A^2) * L/D * abs(Q) puis je crée une matrice diagonale (les K(i) formant donc une diagonale avec les autres éléments de la matrice égaux à 0)

6- On peut alors obtenir une nouvelle matrice M=E' * inv(K) * E qui présente l'avantage de vérifier l'équation suivante:
M*H=C où C désigne la somme des débits à chaque embranchement

7- Il faut alors distinguer les croisements où on connait H (II et III) et ceux où on connaît C (tous les autres).
On va construire deux nouvelles matrices. Pour chaque embranchement i :
A: - Quand C est connu, la colonne i de A est égale à la colonne i de M
- Quand H est connu, l'élément A(i,i) vaut -1 et le reste de la colonne vaut 0
b:Il s'agit un vecteur. Chaque b(i) prend initialement la valeur -1*Somme(M(i,j)*H(j)) (somme des M(i,j)*H pour tous les H connus)
- Quand C est connu, b(i)=b(i)+C
- Quand H est connu, passer à l'élément suivant
NB: J'utilise un vecteur logique (y) pour obtenir chacune de ces deux matrices avec 1 seule ligne de code

8- On peut alors résoudre A*x=b où x est un vecteur contenant pour chaque embranchement:
C si H était connu
H si C était connu

9- Je ré-obtiens alors le vecteur Q via l'équation Q=inv(K)*E*H

10- Je calcule la déviation relative par rapport à l'itération précédente ainsi que le résidu et compare à une tolérance, l'algorithme boucle alors jusqu'à ce que l'une de ces valeurs soit reconnue comme satisfaisante

[wiwaxia]

Bonjour,

Si j'ai bien compris tes explications (sans en avoir saisi tous les détails !), tu recherches par un procédé itératif les limites des 10 composantes du vecteur débit (Q).

Que la convergence soit lente n'est pas étonnant, en raison du grand nombre des variables concernées, et de la non-linéarité des équations.

Peut-être pourrait-tu t'intéresser à la convergence séparée des diverses composantes (q₁, q₂, ... q₁₀): la stabilisation n'est-elle pas plus rapide pour certaines d'entre elles ?
Observe-t-on pour leurs valeurs successives, et à partir d'un certain stade, des suites monotones ou alternées ?

Sur le schéma apparaissent cinq termes dont dépendent tous les autres: (q₁, q₄) et d'autre part (q₂, q₅, q₇), dont la somme correspond à ce qui est fourni par les deux réservoirs:
q₁ + q₄ + q₂ + q₅ + q₇ = 0.2 + 0.3 + 0.5 = 1 m³/s .

Leurs variations ne présenteraient-elles pas, à partir d'un certain nombre d'étapes (k₀) un comportement approximativement géométrique ? Le rapport r_{j, k} = (q_{j, k+2} - q_{j, k+1}) / (q_{j, k+1} - q_{j, k})
deviendrait quasi-constant, et le remplacement toutes les deux étapes du dernier terme obtenu par une combinaison linéaire appropriée:
q'_{j, k+2} = (Ca * q_{j, k+1} + Cb * q_{j, k+2}) permettrait une accélération de la convergence.

Les valeurs: Ca = -r_j/(1 - r_j) , Cb = 1/(1 - r_j) conduisent en théorie directement à la limite;
il sera prudent de retenir pour le rapport une valeur plus modérée, donc plus proche de zéro:
# dans le cas d'une suite monotone (r_j >~ 0.9), prendre r_j = 0.5 , d'où: Ca = -1 , Cb = 2 ;
# dans le cas favorable d'une suite alternée (r_j <~ -0.9), prendre r_j = -0.5 , d'où: Ca = 1/3 , Cb = 2/3 ;
le choix encore plus simple: Ca = Cb = 0.5 (correspondant à r_j = -1) conduit souvent au résultat attendu.

[wiwaxia]

Autre possibilité: si les cinq termes précités résultent de calculs indépendants, leur somme:
S = q₁ + q₄ + q₂ + q₅ + q₇
ne coïncide pas avec la valeur attendue S_tot = 1 m³/s .
Tu pourrais insérer après chaque itération la correction modératrice: q'_j = q_j * (S_tot/S) pour les termes concernés, avant d'en déduire tous les autres.
Mais peut-être ton algorithme en tient-il compte ?

Un procédé analogue est envisageable pour les relations de conservation du débit:
q₅ + q₆ = q₉ (noeud V),
q₇ + q₈ = q₁₀ (noeud VI),
qui consisterait à multiplier les 2 termes de la somme par le facteur f = (q_n / S_lm)^1/2, et celui de droite par son inverse f' = f^-1 .

Il faut exclure la présence de termes négatifs, au risque sinon de ralentir la convergence ou de la faire disparaître; ce qui précède ne saurait s'appliquer aux relations
q₁ + q₂ - q₃ = 0.2 , ou:
q₅ = q₉ - q₆ , par exemple.

Pour dire vrai, je doute un peu de l'efficacité de ces corrections, n'étant pas familier des équations en cause.

[Mani8]

Merci beaucoup pour vos réponses. Je vais regarder tout ça et voir si j'arrive à accélérer un peu le processus mais il est bien possible que je doive me contenter de ce que j'ai déjà.

J'ai cependant un autre problème, le code que je vous ai donné est destiné à résoudre une version simplifiée du système (avec 3 valves fermées). Je dois donc maintenant effectuer la même opération pour le cas de figure où les valves sont ouvertes.
J'ai donc rajouté à mon code les caractéristiques des tuyaux additionnels et ai mis à jour ma matrice de connectivité (celle qui vaut 1 là où un tuyau démarre, -1 là où il finit et 0 ailleurs). Je dois aussi dans ce cas prendre en compte les pertes dues à la résistance des valves (je connais leur valeur K que je me contente donc d'ajouter à celles des tuyaux correspondant).

Cependant, dans ce cas de figure, la formule de Colebrook retourne NaN pour les deux premiers lambda à la 18ème itération.
Je reçois le message d'erreur suivant:
Exiting fzero: aborting search for an interval containing a sign change because complex function value encountered during search.
(Function value at -0.00670718 is -3.51242-13.5481i.)
Check function or try again with a different starting value.


Je ne comprends pas comment cela se fait. La formule fait intervenir la racine du précédent lambda mais j'ai vérifié et celui-ci est bien positif... Je ne vois juste pas d'où peut provenir cette valeur complexe.

EDIT: De toute évidence le problème vient des inputs de fzero, la valeur -0.0067070718 ressemble étrangement à Q(1) pour cette itération c.a.d. 6.8688*10^-4

Le système avec valves ouvertes:


Le code:

Code matlab :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
clear
clear all
close all
 
% Fluid and gravity
rho = 1e3;
mu  = 1e-3;
g   = 9.81;
 
% Colebrook formula
colebrook = @(Re,rr,la) 1/sqrt(la)+2*log10(rr/3.71+2.51/Re/sqrt(la));
 
% Deposit
h1 = 200;
h2 = 200;
 
% Pipe data
L1 = 600;
L2 = 400;
L3 = 900;
L4 = 400;
L6 = 350;
L7 = 600;
L9 = 300;
L10 = 500;
L11 = 300;
L12 = 350;
 
D1 = 0.5;
D2 = 0.5;
D3 = 0.5;
D4 = 0.5;
D6 = 0.5;
D7 = 0.5;
D9 = 0.3;
D10 = 0.3;
D11 = 0.3;
D12 = 0.3;
 
v_L  = [L1 L2 L3 L4 L6 L7 L9 L10 L11 L12];
v_D  = [D1 D2 D3 D4 D6 D7 D9 D10 D11 D12];
 
v_r = 1e-3*ones(1,size(v_L,2));
v_rr = v_r./v_D;
 
% Initial lambda values. Starting with Reynolds infinity.
v_la = 0.25./(log10(v_rr/3.71)).^2;
 
% Initial conditions
v_Q = 0.1*ones(1,size(v_L,2));
v_H = [10 h1 h2 10 10 10 10];
 
%Connectivity matrix
    % 1  2  3  4  6  7  8 
E = [-1  1  0  0  0  0  0; %1
     -1  0  1  0  0  0  0; %2
      1  0  0 -1  0  0  0; %3
      0  1  0 -1  0  0  0; %4
      0  0  1  0 -1  0  0; %6
      0  0  0  1 -1  0  0; %7
      0  0  1  0  0 -1  0; %9
      0  0  0  0  1 -1  0; %10
      0  0  0  0  1  0 -1; %11
      0  0  0  0  0  1 -1];%12
 
 %Declare known C values 
 Unknown=0;
 v_C=[-0.2 Unknown Unknown -0.3 0 0 -0.5];
 y = [0 1 1 0 0 0 0]; %logic matrix: 0 when C is known, 1 when H is known
 
 %%....Main Loop....%%
 
 %Iteration limits
 maxiter = 1000;
 tol = 1e-9;
 inc = 1;
 res = 1;
 iter = 0;
 
 %Creation of solution storage matrices
 Q_results = zeros(size(v_Q,2),maxiter);
 H_results = zeros(size(v_H,2),maxiter);
 inc_results = zeros(1,maxiter);
 res_results = zeros(1,maxiter);
 v_sol = [v_Q v_H(1,1) v_H(1,4:size(v_H,2))];
 
 X=zeros(size(v_H,2),1);
 
 %...........MAIN LOOP.............%
 while (inc>tol) && (res>tol) && (iter<maxiter)
 
 iter = iter+1;
 
 %Storage of previous iteration results
 Q_results(:,iter)=v_Q';
 H_results(:,iter)=v_H';
 v_sol_old = v_sol;
 
 X_old=X;
 
 %Calculate Reynolds numbers from Q
 v_vel = 4*v_Q/pi./v_D.^2;
 v_Re  = abs(rho*v_vel.*v_D/mu);
 
 %Get the new lambdas
 for i = 1:size(v_Q,2)
    v_la(i) = fzero(@(x)colebrook(v_Re(i),v_rr(i),x),v_la(i));
 end
 
 %Obtain matrix M
 v_F = 8/pi^2/g./v_D.^4.*v_la.*v_L./v_D;
 v_K = v_F.*abs(v_Q);
 K = diag(v_K);
 M = E'/K*E;
 
 %Obtain A and b
 A=M*diag(1-y)-diag(y);
 b=-M*(v_H.*y)'+(v_C.*(1-y))';
 
 %Solve the system
 X=A\b;
 
 %Extract the new H and C values 
 v_H=X';
 v_H(2)=h1;
 v_H(3)=h2;
 v_C(2)=X(2);
 v_C(3)=X(3);
 
 %Re-evaluate Q
 v_Q=(K\E*v_H')';
 
 %Update solution matrix
 v_sol = [v_Q v_H(1,1) v_H(1,4:size(v_H,2))];
 
 %Relative increment of solution
 %inc = norm((v_sol-v_sol_old)./v_sol,inf);
 %inc = norm((v_Q-v_sol_old(1:10))./v_Q,inf);
 inc = norm((X-X_old)./X,inf);
 inc_results(iter) = inc;
 
 %Residual of system
 %res = norm(v_sol-v_sol_old,inf);
 res = norm(v_Q-v_sol_old(1:10),inf);
 res_results(iter) = res;
 
 end
 
 %Plots
 inc_results(iter+1:maxiter) = [];
 res_results(iter+1:maxiter) = [];
 
 figure(1)
 subplot(211)
 semilogy(1:iter,inc_results)
 box on,grid on
 xlabel('iteration')
 ylabel('relative increment')
 
 subplot(212)
 semilogy(1:iter,res_results)
 box on,grid on
 xlabel('iteration')
 ylabel('residue')

[wiwaxia]

Il faut s'attendre à ce que le 3^me réservoir, situé nettement plus haut, bouleverse la répartition des flux à l'ouverture des vannes, et inverse le sens de circulation du fluide dans quelques unes des conduites.

Est-il possible de suivre en parallèle les valeurs successives des 14 composantes du vecteur (Q) ? Cela te permettrait de repérer où se produit la première inversion de signe ...
La notation indiquée dans le second schéma est choquante: tu as sans doute de bonnes raisons de l'avoir adoptée, mais j'aurais pour ma part réservé les quatre derniers indices (11, 12 ... 14) aux conduites munies d'une vanne, et dont on peut commander la mise en service.

Et à la réflexion, ce changement de notation est inéluctable, si tu veux maîtriser algorithmiquement l'évolution du réseau d'un état à l'autre (vannes toutes ouvertes ou fermées), avec à la clef l'éventualité de n'en ouvrir que quelques unes ... : les conduites présentes doivent impérativement recevoir les mêmes indices.

Les inversions de signe deviennent à priori gérables au prix d'un petit changement, dans l'expression de la formule de Darcy:

, en remplaçant V² par le produit V*Abs(V) .

PS: le plantage de ton programme provient de la formule de Colebrook, au niveau du nombre de Reynolds: il faut là encore y remplacer (V) par Abs(V).

[Mani8]

Merci beaucoup pour votre aide !!!
J'ai utilisé cette notation pour la simple raison que c'est celle employée dans mon énoncé mais votre remarque est très sensée. Je pense que je vais procéder comme vous le dites, quitte à remettre les résultats dans l'ordre en post-process.