Equation and devices calibration (gates and weirs)

We provide below the scripts of 2 Matlab .m functions (Qouvrage.m and calage_f.m) and a script file .m (calage_r.m) usefull for calibrating a discharge coefficient from mesurement datas using the equations programmed in SIC.

Devices modelisation (gates and weirs) (Script to copy in a file 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

Intermediate function used for calibration (script to copy in a f.m file)

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;

Calibration of discharge coefficient (script to copy in a file 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)

Once these codes have been copied in the .m file as above, you have to run the .m script (calage_r.m) for calculating the discharge coefficient Cd from measures h1, h2, Qm and the description of the device (h1, hv, w, L) using the equations of the evices programmed in SIC (several choices are proposed).