среда, 8 марта 2017 г.

Расход газа при критическом и докритическом истечении (Продувка оборудования - I)


Эта заметка будет слегка нарушать логику изложения – в приведенном в конце заметки листинге будет использована целая куча формул и алгоритмов, ранее в этом блоге не рассматриваемых. Проблема в том, что многие задачи не имеют простого линейного решения. Расчет расхода газа при продувке оборудования входит в число таких задач.
В общем, листинг расчета размещаю без купюр, а вот пояснять его придется в нескольких последующих заметках, разбив пояснения на удобоваримые «куски».
Здесь – про расход газа при критическом и докритическом режимах истечения в атмосферу (формальные признаки режимов истечения были рассмотрены в предыдущей заметке).
Итак.
Вы подходите к свечному крану - пылеуловителя, фильтра-сепаратора, конденсатосборника - и открываете его. Начинается стравливание газа в атмосферу. Эмиссия природного газа, одним словом. Условимся, что давление в источнике газа при этом не снижается, то есть проводится штатная продувка, к примеру, пылеуловителя, находящегося в рабочем режиме.
Чтобы не вводить Вас в заблуждение, будем считать, что длина свечной линии не превышает 10-15 метров, иначе придется учитывать гидравлическое сопротивление трубы, но до задачи расчета пропускной способности трубы мы еще не дошли.
В общем виде расход газа определяется выражением:
Дальше начинаются тонкости.
Вариант первый. Рассчитанная Вами скорость истечения газа из свечной линии меньше, чем скорость звука в газовой среде. Расчет следует выполнять по формуле для докритического режима истечения газа:
Стоит отметить, что режим докритического истечения газа имеет место при перепаде давлений между источником и приемником газа не более (примерно, зависит от давления, температуры и плотности газа) 1,5 кгс/см2.
Вариант второй. Скорость истечения газа равна или превышает скорость звука в газовой среде. В этом случае расчет следует выполнять с учетом критического истечения газа:
Обе зависимости используют одно и то же выражения для расчета коэффициента расхода:
Честно признаюсь – сходу рассчитать расход газа не получится. Как Вы уже обратили внимание, в выражение коэффициента расхода газа (очень эмпирический коэффициент, учитывающий, в том числе, влияние сжатия истекающей струи газа) входит новая характеристика – число Рейнольдса. Его расчет требует отдельного описания, и это будет одной из тем следующих заметок. Как, впрочем, и влияние пропускной способности свечной линии на объемы эмиссии газа. Ну, и количество поворотов 90 градусов на свечной трубе. Промолчу пока про динамическую вязкость, теплопроводность и коэффициент теплопередачи.
Как видите, расчет расхода газа при продувке оборудования и вправду невозможно реализовать простейшими линейными алгоритмами.
Выполнить расчет и «поиграть» величинами давлений, температуры и длиной свечной линии можно в онлайн-режиме на странице сайта:


Серверная часть расчета (PHP):

<?php
// Внимание!
// проверки на адекватность исходных данных
// выполняются на стороне клиента!
//
header('Content-Type^ text/html; charset=utf-8');
if ($_SERVER['HTTP_X_REQUESTED_WITH']=='XMLHttpRequest'){
        if ($_POST){
               // исходные данные из приложения
               $ro = $_POST['ro'];// плотность газа абсолютная, кг/м3
               $pn  = $_POST['pn'];// избыточное давление газа, кгс/см2
               $prt  = $_POST['prt'];// атмосферное давление, мм рт.ст.
               $tn  = $_POST['tn'];// температура газа, по Цельсию
               $tgr  = $_POST['tgr'];// температура грунта, по Цельсию
               $azot  = $_POST['azot'];// молярная составляющая азота
               $ad  = $_POST['ad'];// диаметр свечной линии, мм
               $al  = $_POST['al'];// длина свечной линии, м
               $adk  = $_POST['adk'];// диаметр свечного крана, м
               $bendcount  = $_POST['bendcount'];// количество поворотов 90 градусов
               $strdzeta  = $_POST['strdzeta'];// качество внутренних стенок трубы
               $radius  = $_POST['radius'];// отношение диаметра трубы к радиусу поворота
               $tim  = $_POST['tim'];// время стравливания газа, сек
               // децимальный разделитель
               $ro = str_replace(",", ".", $ro);
               $pn = str_replace(",", ".", $pn);
               $prt = str_replace(",", ".", $prt);
               $tn = str_replace(",", ".", $tn);
               $tgr = str_replace(",", ".", $tgr);
               $azot = str_replace(",", ".", $azot);
               $ad = str_replace(",", ".", $ad);
               $al = str_replace(",", ".", $al);
               $adk = str_replace(",", ".", $adk);
               $bendcount = str_replace(",", ".", $bendcount);
               $radius = str_replace(",", ".", $radius);
               $tim = str_replace(",", ".", $tim);
               // приведение размерности
               $tn = $tn+273.15;// по Кельвину
               $tgr = $tgr+273.15;// по Кельвину
               $azot = $azot/100;// молярная составляющая в долях единицы
               $ad = $ad/1000;// диаметр трубы в метрах
               $adk = $adk/1000;// диаметр свечного крана в метрах
               $patm = $prt*0.001359511;// атмосферное давление в кгс/см2 - для расчета абсолютного давления
               $pk = 0 ;// стравливание в атмосферу - в этом случае
               $tk = $tn-10;// конечная температура
               if (!$al>0) $al = 0.001;// чистое истечение - без учета длины свечной линии
               // загрузка функций
               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_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;    
               }
               // расчет
               // площадь сечения - по меньшему диаметру (труба или кран)
               if ($ad<$adk) {
                       $s = (3.141593*$ad*$ad)/4;
                       $d = $ad;
               } else {
                       $s = (3.141593*$adk*$adk)/4;
                       $d = $adk;
               }
               // промежуточные характеристики
               $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){ // некритическое истечение
                       $pnmpa = ($pn+$patm)*0.0980665;
                       $pkmpa = ($patm)*0.0980665;// конечное давление - атмосферное!
                       $pnpa = $pnmpa*1e6;// Па
                       $pkpa = $pkmpa*1e6;// Па
                       //dmm:=d*1000;// диаметр в мм
                       //
                       $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);
                       //
                       $beta = $pkpa/$pnpa;
                       $q = $kq*$s/$ro*sqrt(2*$k/($k-1)*$pnpa*$ud*(pow($beta, 2/$k)-pow($beta, ($k+1)/$k)))*1;// за 1 сек
                       $qsut = ($q/1000000)*24*60*60;
                       $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);
                       $qexp = $kq*$s/$ro*sqrt(2*$k/($k-1)*$pnpa*$ud*(pow($beta, 2/$k)-pow($beta, ($k+1)/$k)))*$tim;
               } 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)))))*$tim;// м3
               }
               // расчет пропускной способности трубы
               // первое приближение средней температуры
        $km = 1; $cp = 0.65;// начальное предположение
        // предполагаемый транспорт газа - истечение
        $qsut = $qexp/1000000/$tim*60*60*24;// млн.м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));
        // если имеются повороты 90 градусов
        if ($bendcount > 0){
                       if ($strdzeta == 2){
                               $a0 = 0.173022774327032;
                       $a1 = 0.475294874208228;
                       $a2 = -0.833474496518267;
                       $a3 = 0.77102547211251;
                       } else{
                               $a0 = 0.136149068322988;
                       $a1 = 0.043347096743814;
                       $a2 = -0.200992847731959;
                       $a3 = 0.314661561264817;
                       }
                       $dzeta = $a0+$a1*$radius+$a2*$radius*$radius+$a3*$radius*$radius*$radius;
               $leq = $bendcount*($dzeta*$ad/$lambda);// дополнительная эквивалентная длина, м
               $al = $al+$leq;
               //
               $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));
               }
        //
        $q = $qsut*1000000/24/60/60*$tim;// производительность трубы, м3 за указанное время
               if ($q>$qexp){
                       $qres = $qexp;
                       if ($w<$u){$resrt = '<p>Режим истечения газа: некритический</p>';
                       }
                       else {$resrt = '<p>Режим истечения газа: критический</p>';
                       }
               }       else {
                       $qres = $q;
                       $resrt = '<p>Режим истечения газа: ограничен пропускной способностью свечной линии</p>';
               }       
               if ($qres<100){
                       $res = number_format($qres, 3);
                       $res = str_replace(",","",$res);
                       $res = str_replace(".",",",$res);
                       print $resrt.'<p>Расход газа при продувке: '.$res.' м<sup>3</sup></p>';    
               } else {
                       $qres = $qres/1000;
                       $res = number_format($qres, 3);
                       $res = str_replace(",","",$res);
                       $res = str_replace(".",",",$res);
                       print $resrt.'<p>Расход газа при продувке: '.$res.' тыс.м<sup>3</sup></p>';
               }
        }
}
?>