Two-phase filtration in a pipe filled with a porous material.pdf ( -pGg )) = qG; (1) dpGm (1 - 5) IKkG - - div dt H v Математическая модель двухфазной фильтрации в пористых пластах сформулирована на основе анализа работ [1-4]: G..............pG< G G/G\ G 0 P 0 p =p (P ) = p ' —г , P - атмосферное давление; p0 • Уп_ У s - эффективная влагонасыщенность (в принципе меняется от 0 до 1); IK - тензор проницаемости, в общем случае полно заполненный; kG = V1 - s (1 - s2) - относительная проницаемость для газовой фазы; - вязкость газа; qG - источниковый член в балансе газовой фазы; q -источниковый член в балансе жидкости; Pf pF = pF 0 —— , pF 0 - плотность флюида при атмосферном давлении; P kF =yfs (1 - V1 - s2) - относительная проницаемость для жидкой фазы; - вязкость флюида (цF » H° ); ■F— (iKkF dp ms —--div dt G ((F -p)l = qF ; (2) pF' (G \ v p ) - пористость (—, Уп - объем пор); s d - капиллярное давление; p - пороговое pF = pG - pc (s); pc (s) = pl s F d P g давление, причем p =-, где a - параметр, характеризующий пористую среа ду. Таким образом независимыми переменными будут давление газа и влагона-сыщенность. В качестве граничных условий для них будем использовать условия Дирихле. В начальный момент времени давление газа и влагонасыщенность будем считать постоянными и равными значениям на выходе. Одномерное приближение Предельно упростим ситуацию. Примем, что источники отсутствуют, то есть g f " " д q = q = 0. Будем рассматривать одномерное движение, в этом случае V = i —, дх div = —. Кроме того, примем, что движение происходит в горизонтальном надх правлении: (g )х = 0 . В этом случае уравнения (1), (2) можно переписать следующим образом: dpGm (1 - s) д dt дх fдpG ^ KG (3) (4) v дх У/ v Ц KkF dpFms = д F дp дt дх дх V Здесь K - скалярный коэффициент проницаемости. Плотность флюида будем считать неизменной, а плотность газа - линейной функцией давления: G,0 PG =P~0~ pG. В качестве независимых переменных будем рассматривать давление газа и влагонасыщенность. Поэтому давление флюида будем определять по формуле pF = pG - pc (s), где pc (s) = pL 41-s1 . Эта зависимость представлена на рис. 1. s В результате получим два уравнения следующего вида: дpGm (1 - s) = д f KkG (s) G дp (5) дt дх (s))' дх дх дms д f KkF (s) д(p° -p" дt дх Метод контрольного объема Для построения дискретных аналогов исходных дифференциальных уравнений разобьем область решения на контрольные объемы вдоль оси x и проинтегрируем эти уравнения, заменяя объемные интегралы поверхностными. Поскольку стенки трубы непроницаемы, при интегрировании по поверхности рассматриваются только две ее составляющие - два круга, через которые осуществляется течение. Рассматривая значения давления газа и влагонасыщенности в центрах контрольных объемов, а также учитывая все предполагаемые выше упрощения исходных дифференциальных уравнений для случая одномерного приближения, будем иметь разностные уравнения следующего вида: Vim ( (1 - *}) м- ( (1 - *}) - ЙГ • 0ЙГ1 - PSM1 * 0&+1, (7) п+1 - п Vm^Mr- - o^1 -0SS+1- (8) G,n+1 „G,n+1 V С Ъ-° (Z>n+1 \ „G,n+1 K • Sc • kr (sJ ) Pj - Pi где 0;G"+1 --c 'GXJ --, (j - left, right): |aG Ml ->F ,n+1 K • Sc • kF (Sf) pf-1 - pc (s^1)-pG,n+l + pc (sr1) ^ Ml ' Vi - объем i-го контрольного объема, £с - площадь сечений трубы, через которые осуществляется течение, pG n+1 - (pG n1 + pG,n+1)/2 - осредненное значение давления газа на j-й поверхности контрольного объема, pt ,п+ - давление газа в центре /-го контрольного объема, sf+1 - значение насыщенности на j-й поверхности, взятое против потока. Для решения полученной системы 2N дискретных уравнений (N - число контрольных объемов) воспользуемся технологией Ньютона, используемой для решения систем нелинейных уравнений. Выпишем нелинейные невязки уравнений для l-го приближения к величине, изменяемой на временном шаге "+1 в i-м контрольном объеме: (9) right Rk, s vm [ pf,< (i - J)-pf," (i - sf)]-At Z PJJj ; j=left (10) right Rf,i = [s] - sf ] - At Z j . j=left (11) s = 0,3 G рвых s = 0,8 Ркх = 1,2 -105 Па 105 Па Рис. 2. Одномерная область двухфазной фильтрации Определяя производные от невязок по значениям насыщенности и давления газа, запишем для каждого контрольного объема систему линейных алгебраических уравнений из двух уравнений для двух неизвестных приращений рассматриваемых переменных: ( dRf,/ dRf,/1 dp' ds] 5p/ f-rf,/ 1 R,/ l8s] j I- rF ,/ J dp' ds'j J Благодаря тому, что для каждой ячейки записывается независимая система уравнений, легко решаемая с помощью правила Крамера, вычислительный алгоритм становится достаточно простым. Недостатком данного метода являются некоторые ограничения для шага по времени. Вычисления проводились при следующих значениях определяющих параметров: L = 6 м - длина трубы; pF = 1000 кг/м3; р°,° = 1,27 кг/м3; K = 9,869233 • 10-13 м2; = 8,9 • 10-4 Пас; = 1,78 • 10-5 Пас. Далее будем рассматривать движение газо-жидкостной среды в трубе, заполненной пористым материалом, например песком (m = 0,4). Тестирование алгоритма при постоянной насыщенности В случае постоянной насыщенности уравнение для давления газа принимает вид квазилинейного уравнения теплопроводности: Kk_ (s) д ( др_ ~~дГ dp (12) _ дх y_m (1 - s) дх На рис. 3 представлены результаты решения этого уравнения рассматривае мым численным методом. В случае установившегося течения, когда др_ д/ = 0, ре шение легко находится аналитически. На графике аналитическое решение практически совпадает с численным, что говорит о корректности применения рассматриваемого численного метода для решения данных уравнений. + + + аналитическое решение численное решение 36 000 c 108 104 100 _ кПа Рис. 3. Распределения давления газа по длине трубы в моменты времени 36, 108, 180, 360 и 36 000 с В рассматриваемом примере профиль давления газа, изначально имеющий вид «левого нижнего уголка», постепенно деформируется в почти линейное распределение, что характерно для решения уравнения теплопроводности. Численный эксперимент Для анализа аппроксимационной сходимости были проведены вычисления с различными параметрами сетки, а именно шага вдоль оси х. На рис. 4 можно наблюдать характер сходимости решения разностной задачи. Заметной особенностью более точного решения является большая кривизна. В точках перегиба решения совпадают. Рис. 4. Распределение насыщенности в моменты времени 10, 100, 360 часов для различных шагов вдоль оси x (1 - 0,5 м, 2 - 0,2 м, 3 - 0,1 м) При указанных выше параметрах проводился численный эксперимент, отслеживающий динамику развития течения в трубе, заполненной пористым материалом. Шаг сетки вдоль оси x равнялся 0.1 м. В этом случае рассматриваемая область решения была покрыта 60-ю контрольными объёмами. Были получены распределения давления флюида и газа, а также распределение насыщенности флюидом пористого пространства с течением времени, представленное на рис. 5. Полученные результаты позволяют анализировать течение и характер заполнения пористого пространства. Рис. 5. Распределение насыщенности по объему пористого пространства в моменты времени 2, 10, 30, 60, 100, 150, 220, 360 ч Как видно из рис. 5, наблюдается волновой механизм распространения насыщенности по пористому пространству трубы. Причем зона изменения насыщенности с течением времени увеличивается. Характерная выпукло-вогнутая форма ее распределения со временем становится более простой и к моменту, определяющему установившееся состояние течения, кривая s = s(x) имеет вполне определенный знак кривизны. В таблице приводятся данные итерационной сходимости при различных шагах по времени. Точность расчётов - 10-2. Для давления газа столь высокая точность является излишней и приводит к значительному увеличению итераций. Однако эта точность необходима для значений влагонасыщенности, которые колеблятся в пределах от 0 до 1. Данные итерационной сходимости при расчёте течения с общим временем 15 дней Шаг по времени, мин (число шагов) Общее число итераций Среднее число итераций на шаг Время выполнения программы, с 6 (3600) 590 534 768 315 60 (360) 276 477 164 142

Скачать электронную версию публикации