Equation et calage des ouvrages (vannes et seuils)

Nous fournissons ci-dessous les codes de 2 fonctions MatLab .m (Qouvrage.m et calage_f.m) et d’un fichier script .m (calage_r.m) utiles pour caler un coefficient de débit à partir de mesures en utilisant les équations d’ouvrages programmées dans SIC.

Modélisation des ouvrages (vannes et seuils) (code à copier dans un fichier Qouvrage.m)

function [Q,C,w] = Qouvrage(T,h1,h2,w,hv,L,Cd)
% T=1 : Seuil Orifice avec surverse éventuelle (SIC)
%   2 : Seuil Vanne avec surverse éventuelle (SIC)
%   3 : Classique Seuil
%   4 : Classique Orifice
%   5 : Classique Noyee
%   6 : Cunge avec differentes conditions d'ecoulement
% h1  : cote amont (par rapport au radier)
% h2  : cote aval (par rapport au radier)
% w   : ouverture (par rapport au radier)
% hv  : hauteur de la vanne (pour le calcul de la surverse éventuelle)
% L   : largeur
% Cd  : coefficient de debit
% Q   : debit
% C   : code couleur correspondant aux conditions d'ecoulement
%       C=0; debit nul
%       C=1; surface libre denoye
%       C=2; surface libre noye
%       C=3; charge denoye
%       C=4; charge noye partiel
%       C=5; charge noye total
%       C=C+10; si surverse
% Initialisation ---------------------------------------------
R2G = sqrt(2*9.81);
R32 = 3*sqrt(3)/2;
% Tests ------------------------------------------------------
Q=0;
C=0;
if (w==0 & hv==0)
  Q=0;
  return
end
if (h2>h1)
  Q=0;
  return
end
if (h1==0)
  Q=0;
  return
end
% Seuil - Orifice ===========================================
if (T==1 & w>0)
% Conditions d'ecoulement ------------------------------------
  if (h1<=w)
     surfacelibre='oui';
  else
     surfacelibre='non';
  end
  alfa=2/3;
  if (h2<=alfa*h1)
     denoye='oui';
  else
     denoye='non';
     if (h2<=2/3*h1+w/3)
        partiel='oui';
     else
        partiel='non';
     end
  end
% Seuil - Denoye ---------------------------------------------
  if (surfacelibre=='oui') & (denoye=='oui')
     Q=Cd*L*R2G*(h1^1.5);
     C=1;
  end
% Seuil - Noye -----------------------------------------------
  if (surfacelibre=='oui') & (denoye=='non')
     Cs=R32*Cd;
     Q=Cs*L*R2G*((h1-h2)^0.5)*h2;
    C=2;
  end
% Orifice - Denoye -------------------------------------------
  if (surfacelibre=='non') & (denoye=='oui')
     Q=Cd*L*R2G*(h1^1.5-(h1-w)^1.5);
     C=3;
  end
% Orifice - Noye ---------------------------------------------
  if (surfacelibre=='non') & (denoye=='non')
%  Ennoyement partiel ---------------------------------
     if (partiel=='oui')
        Q=Cd*L*R2G*(R32*h2*(h1-h2)^0.5-(h1-w)^1.5);
        C=4;
%  Ennoyement total -----------------------------------
     else
        Cs=R32*Cd;
        Q=Cs*L*R2G*((h1-h2)^0.5)*w;
       C=5;
     end
  end
end
% Seuil - Vanne =============================================
if (T==2 & w>0)
% Calcul de parametres ---------------------------------------
  mu0=2/3*Cd;
% Conditions d'ecoulement ------------------------------------
  if (h1<=w)
     surfacelibre='oui';
     muf=mu0-0.08;
     alfa=0.75;
  else
     surfacelibre='non';
     mu =mu0-0.08/(h1/w);
     mu1=mu0-0.08/(h1/w-1);
     alfa=1-0.14*h2/w;
     if (alfa<0.4)
        alfa=0.4;
     end
     if (alfa>0.75)
        alfa=0.75;
     end
  end
  if (h2<=alfa*h1)
     denoye='oui';
  else
     denoye='non';
     x=sqrt(1-h2/h1);
     beta=-2*alfa+2.6;
     if (x>0.2)
        KF=1-(1-x/sqrt(1-alfa))^beta;
     else
        KF=5*x*(1-(1-0.2/sqrt(1-alfa))^beta);
     end
     alfa1=1-0.14*(h2-w)/w;
     if (alfa1<0.4)
        alfa1=0.4;
     end
     if (alfa1>0.75)
        alfa1=0.75;
     end
     if (h2<=alfa1*h1+(1-alfa1)*w)
        partiel='oui';
     else
        partiel='non';
     end
  end
% Seuil - Denoye ---------------------------------------------
  if (surfacelibre=='oui') & (denoye=='oui')
     Q=muf*L*R2G*(h1^1.5);
     C=1;
  end
% Seuil - Noye -----------------------------------------------
  if (surfacelibre=='oui') & (denoye=='non')
     Q=KF*muf*L*R2G*(h1^1.5);
     C=2;
  end
% Vanne - Denoye ---------------------------------------------
  if (surfacelibre=='non') & (denoye=='oui')
     Q=L*R2G*(mu*h1^1.5-mu1*(h1-w)^1.5);
     C=3;
  end
% Vanne - Noye -----------------------------------------------
  if (surfacelibre=='non') & (denoye=='non')
     x1=sqrt(1-(h2-w)/(h1-w));
     beta1=-2*alfa1+2.6;
     if (x1>0.2)
        KF1=1-(1-x1/sqrt(1-alfa1))^beta1;
     else
        KF1=5*x1*(1-(1-0.2/sqrt(1-alfa1))^beta1);
     end
%  Ennoyement partiel ---------------------------------
     if (partiel=='oui')
        Q=L*R2G*(KF*mu*(h1^1.5)-mu1*(h1-w)^1.5);
        C=4;
%  Ennoyement total -----------------------------------
     else
        Q=L*R2G*(KF*mu*(h1^1.5)-KF1*mu1*(h1-w)^1.5);
        C=5;
     end
  end
end
% Surverse dans cas 1 et 2 ===================================
if hv==0
  hv=inf;
end;
if (T==1) | (T==2)
  if (h1>w+hv)
     surverse='oui';
     alfa=2/3;
     if (h2-w-hv<=alfa*(h1-w-hv))
        denoyesurverse='oui';
     else
        denoyesurverse='non';
     end
  else
     surverse='non';
  end
  % Surverse - Denoye -------------------------------------------
  if (surverse=='oui')
     if (denoyesurverse=='oui')
        Q=Q+0.4*L*R2G*(h1-w-hv)^1.5;
        C=C+10;
     else
        Q=Q+1.04*L*R2G*((h1-h2)^0.5)*(h2-w-hv);
        C=C+10;
     end
  end
end
% Classique - Seuil ==========================================
if (T==3)
  Q=Cd*L*R2G*h1^1.5;
  C=1;
end
% Classique - Orifice ========================================
if (T==4)
  if (h1>w)
     Q=Cd*w*L*R2G*(h1-w/2)^0.5;
     C=3;
  else
     Q=0;
     C=0;
  end
end
% Classique - Noyee ==========================================
if (T==5)
  if (h1>h2)
     Q=Cd*w*L*R2G*(h1-h2)^0.5;
     C=5;
 else
     Q=0;
     C=0;
  end
end
% Classique Cunge ===========================================
if (T==6)
% Conditions d'ecoulement ------------------------------------
  if (h2<=2/3*h1)
     denoye='oui';
     if (w<=2/3*h1)
        surfacelibre='non';
        Q=Cd*L*R2G*w*((h1-w)^0.5);
        C=3;
     else
        surfacelibre='oui';
        Q=Cd*L*R2G/R32*h1^1.5;
        C=1;
     end
  else
     denoye='non';
     if (w<=h2)
        surfacelibre='non';
        Q=Cd*L*R2G*w*((h1-h2)^0.5);
        C=5;
     else
        surfacelibre='oui';
        Q=Cd*L*R2G*h2*((h1-h2)^0.5);
        C=2;
     end
  end
end

Fonction intermédiaire utilisée pour le calage (code à copier dans un f.m)

function Q = f(Cd)
% fonction utilisee pour le calage dans calage\_r.m
% T=1 : Seuil Orifice (SIC)
%   2 : Seuil Vanne (SIC)
%   3 : Classique Seuil
%   4 : Classique Orifice
%   5 : Classique Noyee
% h1  : cote amont (par rapport au radier)
% h2  : cote aval (par rapport au radier)
% w   : ouverture
% hv  : hauteur de la vanne (pour le calcul de la surverse éventuelle)
% L   : largeur
% Cd  : coefficient de debit
% Q   : debit calcule par Qouvrage
% Qm  : debit mesure pour le calage
% C   : code de couleur correspondant aux conditions d'ecoulement
% Initialisation
T = 2;
hr = 41.95
h1 = 43.28 - hr;
h2 = 42.77 - hr;
hv = 0;
w = 0.66;
L = 1.39;
Qm = 1.909;
% Calcul
[Q,C,w] = Qouvrage(T,h1,h2,w,hv,L,Cd);
Q = Q - Qm;

Calage du coefficient de débit (code à copier dans un fichier calage_r.m)

% Fichier .M utilisé pour caler Cd
% fzero est une fonction MatLab qui trouve le zéro
% d'une fonction proche d'une valeur donnée
Cd = fzero('calage_f',0.5)

Une fois que ces codes ont été copiés dans les fichiers .m comme indiqué ci-dessus, il suffit d’exécuter le fichier script .m (calage_r.m) pour calculer le coefficient de débit Cd à partir de mesures h1, h2, Qm et de la description de l’ouvrage (hr, hv, w, L) en utilisant les équations d’ouvrages programmées dans SIC (divers choix sont proposés).