Changeset 60828 in spip-zone


Ignore:
Timestamp:
Apr 30, 2012, 6:45:46 PM (9 years ago)
Author:
dorch@…
Message:

Développement des calculs d'ouvrages. Toutes les paramètres peuvent être en calcul. Lois 1,3,4,5 programmées.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • _plugins_/hydraulic/branches/v0.3/hyd_inc/ouvrage.class.php

    r60450 r60828  
    5353     * - CT : Coefficient de débit partie triangulaire pour les trapézoïdales
    5454     * - CS : Coefficient de débit de la surverse
     55     * - P : Précision du calcul
    5556     */
    5657    private $tP = array();
    5758
    58     const G = 9.81;
    59     const R2G = 4.42944; //sqrt(2*self::gP);
    60     const R32 = 2.59807; // 3*sqrt(3)/2;
     59    const G = 9.81; /// Constante de gravité terrestre
     60    const R2G = 4.42944; /// sqrt(2*self::gP);
     61    const R32 = 2.59807; /// 3*sqrt(3)/2;
     62    const IDEFINT = 100; /// Pas de parcours de l'intervalle pour initialisation dichotomie
     63    const IDICMAX = 100; /// Itérations maximum de la dichotomie
    6164
    6265    /**
     
    104107     * @return array(0=> donnée calculée, 1=> Flag d'écoulement)
    105108     */
    106     public function Calc($sCalc,$rInit=0) {
    107         $rVarC = &$this->tP[$sCalc];
    108 
     109    public function Calc($sCalc,$rInit=0.) {
     110        //print_r($this->tP);
    109111        // Calcul du débit (facile !)
    110         if($sCalc=='Q') return $this->OuvrageQ();
    111 
    112         // Sinon calcul d'une autre donnée par dichotomie
    113 
     112        if($sCalc=='Q') {
     113            return $this->OuvrageQ();
     114        }
     115        else {
     116            // Sinon calcul d'une autre donnée par dichotomie
     117            $rVarC = &$this->tP[$sCalc];
     118            $QT = $this->tP['Q']; // Débit recherché (Target)
     119            $XMinInit = 0;
     120            $rVarC = $XMinInit;
     121            list($Q1,$nFlag) = $this->OuvrageQ();
     122            $XMaxInit = $rInit*10; /// @todo Boucler la valeur max sur 10,100,1000,10000
     123            $rVarC = $XMaxInit;
     124            list($Q2,$nFlag) = $this->OuvrageQ();
     125            $DX = ($XMaxInit - $XMinInit) / floatval(self::IDEFINT);
     126            $nIterMax = floor(max($XMaxInit - $rInit,$rInit - $XMinInit) / $DX + 1);
     127            $Xmin = $rInit;
     128            $Xmax = $rInit;
     129            $X1 = $rInit;
     130            $X2 = $rInit;
     131            $rVarC = $rInit;
     132            list($Q,$nFlag) = $this->OuvrageQ();
     133            $Q1 = $Q;
     134            $Q2 = $Q;
     135            //echo "\n".'nIterMax='.$nIterMax.'  XMinInit='.$XMinInit.'  XMaxInit='.$XMaxInit.'  DX='.$DX;
     136
     137
     138            for($nIter=1;$nIter<=$nIterMax;$nIter++) {
     139                //Ouverture de l'intervalle des deux côtés puis à droite et à gauche
     140                $Xmax = $Xmax + $DX;
     141                if($Xmax > $XMaxInit xor $DX <= 0) $Xmax = $XMaxInit;
     142                $rVarC = $Xmax;
     143                list($Q,$nFlag) = $this->OuvrageQ();
     144                if($Q1 < $Q2 xor $Q <= $Q2) {
     145                    $Q2 = $Q;
     146                    $X2 = $Xmax;
     147                }
     148                if($Q1 < $Q2 xor $Q >= $Q1) {
     149                    $Q1 = $Q;
     150                    $X1 = $Xmax;
     151                }
     152                $Xmin = $Xmin - $DX;
     153                if($Xmin < $XMinInit xor $DX <= 0) {
     154                    $Xmin = $XMinInit;
     155                }
     156                $rVarC = $Xmin;
     157                list($Q,$nFlag) = $this->OuvrageQ();
     158                if($Q1 < $Q2 xor $Q <= $Q2) {
     159                    $Q2 = $Q;
     160                    $X2 = $Xmin;
     161                }
     162                if($Q1 < $Q2 xor $Q >= $Q1) {
     163                    $Q1 = $Q;
     164                    $X1 = $Xmin;
     165                }
     166/*
     167                echo "\n".'nIter='.$nIter.' Xmin='.$Xmin.' Xmax='.$Xmax;
     168                echo "\n".'X1='.$X1.' Q1='.$Q1.' X2='.$X2.' Q2='.$Q2;
     169                echo "\n".'$QT > $Q1 xor $QT >= $Q2 = '.($QT > $Q1 xor $QT >= $Q2);
     170*/
     171                if($QT > $Q1 xor $QT >= $Q2) {break;}
     172            }
     173
     174            if($nIter >= self::IDEFINT) {
     175                // Pas d'intervalle trouvé avec au moins une solution
     176                if($Q2 < $QT and $Q1 < $QT) {
     177                    // Cote de l'eau trop basse pour passer le débit il faut ouvrir un autre ouvrage
     178                    $rVarC = $XmaxInit;
     179                }
     180                else {
     181                    // Cote de l'eau trop grande il faut fermer l'ouvrage
     182                    $rVarC = $XminInit;
     183
     184                }
     185                list($Q,$nFlag) = $this->OuvrageQ();
     186                $nFlag = -1;
     187            }
     188            else {
     189                // Dichotomie
     190                $X = $rInit;
     191                for($nIter = 1; $nIter<=self::IDICMAX;$nIter++) {
     192                    $rVarC=$X;
     193                    list($Q,$nFlag) = $this->OuvrageQ();
     194                    if(abs($Q/$QT-1.) <= $this->tP['P']) {break;}
     195                    if($QT < $Q xor $Q1 <= $Q2) {
     196                        // QT < IQ et Q(X1) > Q(X2) ou pareil en inversant les inégalités
     197                        $X1=$rVarC;
     198                    }
     199                    else {
     200                        // QT < IQ et Q(X1) < Q(X2) ou pareil en inversant les inégalités
     201                        $X2=$rVarC;
     202                    }
     203                    $X=($X2+$X1)*0.5;
     204                }
     205                if($nIter == self::IDICMAX) {
     206                    //IF1 <-- -10 anomalie: la dichotomie n'a pas abouti en ITER iterations
     207                    $nFlag = -1;
     208                }
     209            }
     210        }
    114211        return array($rVarC,$nFlag);
    115212    }
     
    180277            $this->tP['ZM'] = $ZM;
    181278        }
    182 
     279        //echo "\n".'OuvrageQ='.$rQ.' / '.$nFlag;
    183280        return array($rQ,$nFlag);
    184281    }
Note: See TracChangeset for help on using the changeset viewer.