пятница, 10 марта 2017 г.

Пропускная способность участка газопровода и коэффициент гидравлической эффективности


Пропускная способность участка магистрального газопровода – очень важная характеристика как для проектировщиков газотранспортной системы, так и для производственно-диспетчерских служб. Рассчитывается она следующим образом:
Под разницей высот имеется в виду разница между высотой начала и конца участка магистрального газопровода над уровнем моря.
В дополнение к изложенным в предыдущих заметках зависимостям, привожу формулу для расчета коэффициента гидравлического сопротивления:
Все занимательное в этом расчете кроется в алгоритме, который (вот удивительно!) замкнут сам на себя.
Итак, вариант алгоритма для расчета пропускной способности участка магистрального газопровода:
Первый этап:
  • Рассчитать относительную плотность газа;
  • Рассчитать среднее давление и среднюю температуру газа. Среднюю температуру рассчитывать по упрощенной формуле. Если температура газа в конце участка ниже температуры грунта – средняя температура равна температуре газа в конце участка газопровода;
  • Рассчитать коэффициент сжимаемости для средних значений давления и температуры газа;
  • Сделать предположение об объеме транспорта газа. На данном этапе это может быть любая разумная величина, точность пока не важна. В реале у меня прописана функция, привязанная к диаметрам трубы – я описывал ее в алгоритме для расчета числа Рейнолдса.
Второй этап:
  • Рассчитать динамическую вязкость. При расчете приведенных значений давления и температуры использовать средние значения параметров;
  • Рассчитать число Рейнолдса;
  • Рассчитать коэффициент гидравлического сопротивления λ;
  • Рассчитать теоретическую пропускную способность трубы.
Второй этап расчета нужно выполнить повторно, использовав для пересчета числа Рейнолдса  полученное значение пропускной способности трубы.
Получено более-менее приближенное к реальности значение пропускной способности участка магистрального газопровода, поэтому можно перейти к третьему этапу расчета:
  • По точной формуле рассчитать среднюю температуру газа. При расчете использовать последнее из рассчитанных значений пропускной способности трубы;
  • Рассчитать коэффициент сжимаемости газа;
  • Рассчитать число Рейнолдса;
  • Рассчитать коэффициент гидравлического сопротивления λ;
  • Рассчитать теоретическую пропускную способность трубы.
Третий этап расчета необходимо выполнить два раза.
Расчет закончен.
Расчет коэффициента гидравлической эффективности участка газопровода выполняется по следующей формуле:
Таким образом, перед расчетом коэффициента гидравлической эффективности необходимо выполнить предыдущий расчет - пропускной способности трубы, используя в качестве предположительного транспорта газа его фактическое значение.
Отмечу, что расчет пропускной способности трубы (впрочем, как и прочие расчеты, связанные с определением коэффициента теплопередачи) выполняется для подземного исполнения газопровода. Для надземного исполнения, которое, к примеру, имеет место в ОАО «Норильскгазпром», приведенные здесь расчетные методики потребуют значительной коррекции.
Небольшое примечание. Помимо пропускной способности, в результат расчета очень полезно выводить заодно и линейную скорость потока газа в начале и в конце участка газопровода. Очень полезно для контроля допустимых 20 м/сек.
Онлайн-расчет пропускной способности участка магистрального газопровода:


Онлайн-расчет коэффициента гидравлической эффективности:


Исходный код серверной части расчета (PHP) пропускной способности участка газопровода:
<?php
// Внимание!
// проверки на адекватность исходных данных
// выполняются на стороне клиента!
//
header('Content-Type^ text/html; charset=utf-8');
if ($_SERVER['HTTP_X_REQUESTED_WITH']=='XMLHttpRequest'){
        if ($_POST){
               // исходные данные из приложения
               $ro = $_POST['ro'];// плотность газа абсолютная, кг/м3
               $azot  = $_POST['azot'];// молярная составляющая азота
               $pn  = $_POST['pn'];// избыточное давление газа в начале МГ, кгс/см2
               $pk  = $_POST['pk'];// избыточное давление газа в конце МГ, кгс/см2
               $prt  = $_POST['prt'];// атмосферное давление, мм рт.ст.
               $tn  = $_POST['tn'];// температура газа в начале МГ, по Цельсию
               $tk  = $_POST['tk'];// температура газа в конце МГ, по Цельсию
               $tgr  = $_POST['tgr'];// температура грунта, по Цельсию
               $ad  = $_POST['ad'];// диаметр МГ, мм
               $al  = $_POST['al'];// длина МГ, км
               $hh1  = $_POST['hh1'];// начало МГ - высота над уровнем моря, м
               $hh2  = $_POST['hh2'];// конец МГ - высота над уровнем моря, м
               // децимальный разделитель
               $ro = str_replace(",", ".", $ro);
               $azot = str_replace(",", ".", $azot);
               $pn = str_replace(",", ".", $pn);
               $pk = str_replace(",", ".", $pk);
               $prt = str_replace(",", ".", $prt);
               $tn = str_replace(",", ".", $tn);
               $tk = str_replace(",", ".", $tk);
               $tgr = str_replace(",", ".", $tgr);
               $ad = str_replace(",", ".", $ad);
               $al = str_replace(",", ".", $al);
               $hh1 = str_replace(",", ".", $hh1);
               $hh2 = str_replace(",", ".", $hh2);
               // приведение размерности
               $tn = $tn+273.15;// по Кельвину
                $tk = $tk+273.15;// по Кельвину
               $tgr = $tgr+273.15;// по Кельвину
               $azot = $azot/100;// молярная составляющая в долях единицы
               $ad = $ad/1000;// диаметр трубы в метрах
               $al = $al*1000;// длина участка МГ в метрах
               //
               $patm = $prt*0.001359511;// атмосферное давление в кгс/см2
               // загрузка функций
               function f_delta($ro) {
               // относительная плотность газа
               // ro - плотность газа, кг/м3
               $delta = $ro/1.2044;
               return $delta;
               }
               function f_psr($pn, $pk, $patm){
               // среднее избыточное давление на участке МГ
               // pn - начальное избыточное давление, кгс/см2
               // pk - конечное избыточное давление, кгс/см2
               $pna = $pn+$patm;
               $pka = $pk+$patm;
               $psr = 2/3*($pna+($pka*$pka)/($pna+$pka))-$patm;
               return $psr;   
               }
               function f_z($delta, $pn, $patm, $tn){
               // коэффициент сжимаемости газа
               // delta - дельта
               // pn - избыточное давление газа, кгс/см2
               // patm - атмосферное давление, кгс/см2
               // tn - температура газа, по Кельвину
               $pabs = $pn+$patm;// абсолютное давление
               $pmpa = $pabs*0.0980665;// перевод кгс/см2 в МПа
               $z = 1-((10.2*$pmpa-6)*(0.345/100*$delta-0.446/1000)+0.015)*(1.3-0.0144*($tn-283.2));
               return $z;     
               }
               function f_adiabata($azot, $ro, $pn, $patm, $tn){
               // коэффициент адиабаты
               // azot - молярная составляющая азота, доли единицы
               // ro - плотность газа абсолютная, кг/м3
               // pn - избыточное давление, кгс/см2
               // patm - атмосферное давление, кгс/см2
               // tn - температура, по Кельвину
               $pabs = $pn+$patm;// абсолютное давление
               $pmpa = $pabs*0.0980665;// перевод кгс/см2 в МПа
               $k1 = 1.556*(1+0.074*$azot);
               $k2 = 3.9/10000*$tn*(1-0.68*$azot)+0.208*$ro;// последний плюс - не ошибка! потом будем отнимать все сразу
               $k3 = pow(($pmpa/$tn), 1.43)*(384*(1-$azot)*(pow(($pmpa/$tn), 0.8))+26.4*$azot);
               $k = $k1-$k2+$k3;// (Сарданашвили, стр. 44, формула 1.69)
               return $k;     
               }
               function f_rofact($ro, $pn, $patm, $tn, $z){
               // фактическая плотность газа
               // ro - абсолютная плотность газа
               // pn - избыточное давление, кгс/см2
               // patm - атмосферное давление, кгс/см2
               // tn - температура газа, по Кельвину
               // z - коэффициент сжимаемости
               $pabs = $pn+$patm;// абсолютное давление
               $pmpa = $pabs*0.0980665;// перевод кгс/см2 в МПа
               $ud = $ro*(293.15*$pmpa)/($tn*0.101325*$z);
               return $ud;    
               }
               function f_speedsound($tn, $adiabata, $z, $ro){
               // скорость звука в газе, м/с
               // tn - температура газа, по Кельвину
               // adiabata - коэффициент адиабаты
               // z - коэффициент сжимаемости газа
               // ro - плотность газа абсолютная, кг/м3
               $a = ($tn*$adiabata*$z)/$ro;
               $u = 18.591*pow($a, 0.5);
               return $u;     
               }
               function f_speedgas($adiabata, $rofact, $pn, $patm, $pk){
               // скорость истечения газа (по формуле Сен-Венана), м/с
               // adiabata - коэффициент адиабаты
               // rofact - фактическая плотность газа, кг/м3
               // pn - начальное избыточное давление, кгс/см2
               // pk - конечное избыточное давление, кгс/см2
               $pnmpa = ($pn+$patm)*0.0980665;
               $pkmpa = ($pk+$patm)*0.0980665;// внимание! если в атмосферу - то ноль в исходные данные
               $a = $pkmpa/$pnmpa;
               $b =($adiabata-1)/$adiabata; 
               $w = sqrt(2*$adiabata/($adiabata-1)*$pnmpa*1000000/$rofact*(1-pow($a, $b)));
               return $w;     
               }
               function f_qf($d){
               // начальное предположение о пропускной способности газа, млн.м3/сут
               // d - диаметр трубы в метрах
               $ka = 0.138177961111069;
               $kb = -2.48613691329956;
               $kc = 13.8517427444458;
               $kd = 11.1694717407227;
               //
               $qsut = $ka+$kb*$d+$kc*$d*$d+$kd*$d*$d*$d;
               return $qsut;
               }
               function f_mu($ro, $pn, $tn){
               // динамический коэффициент вязкости, кгс*с/м2
               $ppk = 1.808*(26.831-1.205*($ro/1.2044));
               $tpk = 155.24*(0.564+1.205*($ro/1.2044));
               $ppr = $pn/$ppk;
               $tpr = $tn/$tpk;
               $musi = 5.1e-6*(1+$ro*(1.1-0.25*$ro))*(0.037+$tpr*(1-0.104*$tpr))*(1+($ppr*$ppr)/(30*($tpr-1)));
               $mu = $musi/9.81;
               return $mu;    
               }
               function f_re($ro, $qsut, $d, $mu){
               // число Рейнолдса
               // qsut - расход газа, млн.м3/сут
               // ro - плотность газа, кг/м3
               // d -  внутренний диаметр, м
               // mu - динамическая вязкость
               $dmm = $d*1000;
               $re = 1810*$qsut*($ro/1.2044)/($dmm*$mu);
               return $re;    
               }
               function f_expnocrit($s, $tim, $pn, $pk){
               // некритическое истечение
               // s - площадь истечения, м2
               // tim - время истечения, сек
               // pn - начальное давление, кгс/см2
               // pk - конечное давление, кгс/см2
               $dp = $pn-$pk;
               //
               $qexp = 110*$s*$dp*$tim;
               return $qexp;  
               }
               function f_ax($km, $d, $l, $q, $delta, $cp){
               // множитель для определения средней температуры
               // km - коэффициент теплопроводности
               // d - внутренний диаметр, м
               // l - длина МГ, м
               // q - транспорт газа, млн.м3/сут
               // delta - относительная плотность газа
               $dmm = $d*1000;
               $lkm = $l/1000;
               $ax = 62.6*(($km*$dmm*$lkm)/($q*$delta*$cp*1e6));
               return $ax;    
               }
               function f_tsr($ax, $tgr, $tn, $tk){
               // средняя температура газа, по Кельвину
               // ax - множитель для Trs
               // tgr - температура грунта, по Кельвину
               // tn - температура в начале МГ, по Кельвину
               // tk - температура в конце МГ, по Кельвину   
               //$tsr = $tgr+(($tn-$tk)/$ax)*(1-exp(-$ax));
               if ($tgr<$tk){
               $tsr = $tgr+($tn-$tk)/log(($tn-$tgr)/($tk-$tgr));}
               else {$tsr = $tk;}
               return $tsr;
               }
               function f_cp($tsr, $psr){
               // теплоемкость газа в смешанной системе координат, кКал/кг*К
               // tsr - средняя температура, по Кельвину
               // psr - среднее давление, кгс/см2
               $cp = 0.405+4.39e-4*$tsr+46000*($psr-1)/($tsr*$tsr*$tsr);
               return $cp;    
               }
               function f_di($tsr, $cp){
               // коэффициент Джоуля-Томсона, К/(кгс/см2)
               // tsr - средняя температура, по Кельвину
               // cp - теплоемкость газа
               $di = (23000/($tsr*$tsr)-0.035)/$cp;
               return $di;    
               }
               function f_lambda($re, $d){
               // лямбда- коэффициент сопротивления трения
               // re - число Рейнолдса
               // d - внутренний диаметр, м
               $dmm = $d*1000;
               $lambda = 0.067*(exp(log(158/$re+0.06/$dmm)*(1/5)));
               return $lambda; 
               }
               function f_eax($tn, $tk, $tgr, $pn, $pk, $psr, $patm, $ax, $di){
               // множитель для Km
               // tn - начальная температура, по Кельвину
               // tk - конечная температура, по Кельвину
               // tgr - температура грунта, по Кельвину
               // pn - начальное избыточное давление, кгс/см2
               // pk - конечное избыточное давление, кгс/см2
               // psr - среднее избыточное давление, кгс/см2
               // ax - множитель для Tsr
               // di - коэффициент Джоуля-Томсона
               $pnabs = $pn+$patm;
               $pkabs = $pk+$patm;
               $psrabs = $psr+$patm;
               $deltat = abs($tn);
               for ($x=1; $x<=3000; $x++){
                       $dx = ($ax/100)*$x;
                       $tcrush = $tgr+($tn-$tgr)*exp(-$dx)-$di*($pnabs*$pnabs-$pkabs*$pkabs)/(2*$dx*$psrabs)*(1-exp(-$dx));
                       if ($deltat>abs($tcrush-$tk)){
                               $eax = $dx;
                       $deltat = abs($tcrush-$tk);
                       }
               }
               return $eax;   
               }
               function f_kmdross($eax, $q, $delta, $cp, $d, $l){
               // коэффициент теплопередачи с учетом эффекта дросселирования, кКал/м2чК
               // eax - множитель для KM
               // q - транспорт газа, млн.м3/сут
               // delta - относительная плотность газа
               // cp - теплоемкость газа
               // d - диаметр МГ, м
               // l - длина МГ, м
               $dmm = $d*1000;
               $lkm = $l/1000;
               $km = ($eax*$q*$delta*$cp*1e6)/(62.6*$dmm*$lkm);
               return $km;    
               }
               // расчет
               // площадь сечения 
               $s = (3.141593*$ad*$ad)/4;
               $d = $ad;
               // промежуточные характеристики
               $delta = f_delta($ro);
               $z = f_z($delta, $pn, $patm, $tn);// по условиям МГ
               $k = f_adiabata($azot, $ro, $pn, $patm, $tn);
               $ud = f_rofact($ro, $pn, $patm, $tn, $z);
               $u = f_speedsound($tn, $k, $z, $ro);
               $w = f_speedgas($k, $ud, $pn, $patm, $pk);
               // сравнение скорости звука и скорости истечения газа
               if ($w<$u){ // некритическое истечение
                       $qexp = f_expnocrit($s, 1, $pn, $pk);// м3
               } else { // критическое истечение - через массовый расход
                       // перевод давлений в МПа
                       $pnmpa = ($pn+$patm)*0.0980665;
                       $pkmpa = ($patm)*0.0980665;// конечное давление - атмосферное!
                       $pnpa = $pnmpa*1e6;// Па
                       // первая итерация:
                       $qsut = f_qf($d);// предположительный расход газа, млн.м3/сут
                       $mu = f_mu($ro, $pn, $tn);// динамическая вязкость
                       $re = f_re($ro, $qsut, $d, $mu);// число Рейнолдса
                       // коэффициент расхода:
                       $kq = 0.587+(5.5/(pow($re, 0.5)))+(0.348/(pow($re, 0.333)))-(110.92/$re);
                       // объем стравленного газа за 1 секунду:
                       $q = $kq*$s/$ro*sqrt($k*$pnpa*$ud*(pow((2/($k+1)), (($k+1)/($k-1)))))*1;
                       // перевод расхода за секунду в суточный и млн.м3:
                       $qsut = ($q/1000000)*24*60*60;// млн.м3/сут
                       // вторая итерация:
                       $re = f_re($ro, $qsut, $d, $mu);// число Рейнолдса
                       // коэффициент расхода:
                       $kq = 0.587+(5.5/(pow($re, 1/2)))+(0.348/(pow($re, 0.333)))-(110.92/$re);
                       // объем стравленного газа за tim секунд:
                       $qexp = $kq*$s/$ro*sqrt($k*$pnpa*$ud*(pow((2/($k+1)), (($k+1)/($k-1)))))*1;// м3
               }
               // расчет пропускной способности трубы
               // первое приближение средней температуры
        $km = 1; $cp = 0.65;// начальное предположение
        // предполагаемый транспорт газа - истечение
        //$qsut = $qexp/1000000/1*60*60*24;// млн.м3/сут
               $qsut = f_qf($d);// предположительный расход газа, млн.м3/сут
        $ax = f_ax($km, $ad, $al, $qsut, $delta, $cp);// множитель для Tsr
        $tsr = f_tsr($ax, $tgr, $tn, $tk);// средняя температура
               $psr = f_psr($pn, $pk, $patm);
        $cp = f_cp($tsr, $psr);// теплоемкость газа
        $di = f_di($tsr, $cp);// коэффициент Джоуля-Томсона
        $eax = f_eax($tn, $tk, $tgr, $pn, $pk, $psr, $patm, $ax, $di);
        $km = f_kmdross($eax, $qsut, $delta, $cp, $ad, $al);
        // второе приближение средней температуры
        $ax = f_ax($km, $ad, $al, $qsut, $delta, $cp);// множитель для Tsr
        $tsr = f_tsr($ax, $tgr, $tn, $tk);// средняя температура
        // теоретический транспорт газа
        $z =f_z($delta, $psr, $patm, $tsr);// коэф. сжимаемости
        $mu = f_mu($ro, $psr, $tsr);// динамическая вязкость
        $re = f_re($ro, $qsut, $ad, $mu);// число Рейнолдса
        $lambda = f_lambda($re, $ad);// лямбда
               // приведение размерностей для расчета расхода
               $lkm = $al/1000;// длина трубы в км
        $dmm = $ad*1000;// диаметр трубы в мм
               $pnabs = $pn+$patm;
               $pkabs = $pk+$patm;
        $qsut =3.26e-7*exp(log($dmm)*(5/2))*sqrt(($pnabs*$pnabs-$pkabs*$pkabs)/($delta*$lambda*$z*$tsr*$lkm));// теор. расход промежуточный
               // транспорт газа с учетом высоты над уровнем моря
        $re = f_re($ro, $qsut, $ad, $mu);// число Рейнолдса
        $lambda = f_lambda($re, $ad);// лямбда
        $dH = $hh2-$hh1;
        $A = 62.6*$km*$dmm/($qsut*$delta*$cp*1000000);// множительная чепухня
        $AL = $A*$lkm;
        $am = $delta/(14.64*$tsr*$z);
        $b = 1+$am*$dH/2;
        $qsut = 3.26e-7*exp(log($dmm)*(5/2))*sqrt(($pnabs*$pnabs-$pkabs*$pkabs*(1+$am*$dH))/($delta*$lambda*$z*$tsr*$lkm*$b));// теор. расход с учетом dH
               $qtheor = $qsut/24*1000;// теоретический расход - тыс.м3/час
               // граничные условия проверки на критический режим истечения
               $deltap = $pn-$pk;
               if ($deltap>15 AND $lkm<1){// разница давлений и длина ЛЧ
                       $qstrav = $qexp*60*60/1000;// стравливание - тыс.м3/час
               } else {$qstrav = 0;}
               if ($qstrav>0 AND $qstrav<$qtheor){$res = $qstrav;}
               else {$res = $qtheor;}
               //
               $qsut = $res*24/1000;
               $s = 3.14159*$ad*$ad/4;
               // начало МГ
               $z = f_z($delta, $pn, $patm, $tn);// по условиям МГ
               $ud = f_rofact($ro, $pn, $patm, $tn, $z);
               $qmin = 2.45*$qsut*$z*$tn/$pnabs;// расход - м3/мин - фактический, т.е. сжатый газ
               $v = $qmin/$s;// линейная скорость - м/мин
                $vmn = $v/60;// линейная скорость - м/сек
               // конец МГ
               $z = f_z($delta, $pk, $patm, $tk);// по условиям МГ
               $ud = f_rofact($ro, $pk, $patm, $tk, $z);
               $qmin = 2.45*$qsut*$z*$tk/$pkabs;// расход - м3/мин - фактический, т.е. сжатый газ
               $v = $qmin/$s;// линейная скорость - м/мин
               $vmk = $v/60;// линейная скорость - м/сек
        if ($res>0){
                       $res = number_format($res, 3);
                       $res = str_replace(",","",$res);
                       $res = str_replace(".",",",$res);
                       //
                       $km = number_format($km, 3);
                       $km = str_replace(",","",$km);
                       $km = str_replace(".",",",$km);
                       //
                       $vmn = number_format($vmn, 2);
                       $vmn = str_replace(",","",$vmn);
                       $vmn = str_replace(".",",",$vmn);
                       //
                       $vmk = number_format($vmk, 2);
                       $vmk = str_replace(",","",$vmk);
                       $vmk = str_replace(".",",",$vmk);
                       //
                       print '<p>Расчетная пропускная способность: '.$res.' тыс.м<sup>3</sup>/час</p>';
                       print '<p>Предположительный коэффициент теплопередачи: '.$km.' кКал/м<sup>2</sup>чК</p>';
                       print '<p>Линейная скорость потока газа:</p>';
                       print '<p>В начале участка газопровода: '.$vmn.' м/сек</p>';
                       print '<p>В конце участка газопровода: '.$vmk.' м/сек</p>';
               } else {
                       print '<p>Ошибка расчета: нулевое расчетное значение транспорта газа</p>';
               }
        }
}
?>