Введение
Повышение эффективности добычи углеводородов является важной научно-технической задачей нефтегазовой отрасли. В немалой степени это обусловлено истощением традиционных месторождений, а также увеличивающейся потребностью в энергии. Тем не менее ясно, что увеличение добычи не может быть достигнуто без возможности гибкого и своевременного управления месторождением [1, 2], что подразумевает наличие его прогнозной модели, на основе которой и будут приниматься все управленческие решения.
Существуют классические подходы к построению таких моделей, например с использованием известных физических и эмпирических законов, позволяющих представить некоторый физический процесс в виде системы дифференциальных и алгебраических уравнений. С другой стороны, доступность бюджетных датчиков давления и температуры значительным образом увеличила количество данных, получаемых со скважины. В данной работе предлагается гибридный подход к моделированию, в котором общая форма уравнений количества движения определяется на основе нестационарных данных, измеряемых в определенном временном интервале для ограниченного числа пространственных точек. Неопределенные физические параметры модели затем уточняются при помощи методов оптимизации. Хотя предлагаемые в настоящей работе методы носят общий характер, целью исследования является их применение к задаче идентификации динамической модели течения многофазного флюида в скважине.
Метод
Рассмотрим традиционный гидродинамический подход [3] к описанию течения двухфазного флюида в скважине длиной L, наклоненной под углом θ (рис. 1), основанный на использовании одномерных уравнений и полуэмпирических замыкающих соотношений. В частности, используется модель потока дрейфа [4], в рамках которой рассматривается движение смеси фаз как единого целого. Тогда имеем уравнение динамики в виде:
. (1)
Тогда течение флюида характеризуется следующими величинами, зависящими от координаты x вдоль оси скважины и от времени t: скоростями vg газа и vl нефти, долями αg и αl, занимаемыми газом и нефтью в площади сечения трубы, а также давлением p. Под v понимается скорость смеси флюидов, представляющая собой отношение общего объемного расхода смеси к площади поперечного сечения трубы. Приближенно нефть считается несжимаемой, а для газа выполняется идеальный закон. Последний член в правой части уравнения (1) отвечает за трение флюида о стенки трубы. Входящий в него параметр k, вообще говоря, зависит от целого ряда факторов, вплоть до числа Рейнольдса, характеризующего течение. Математическое описание процесса течения дополняется законом сохранения масс для каждой фазы, а также моделью проскальзывания фаз [5], которая связывает скорость смеси со скоростью газа.
. (2)
При этом величины C и vd в общем случае представляют собой довольно сложные функции параметров течения [5].
Рисунок 1. Схематическое изображение скважины (или ее участка) длиной L. Голубым цветом показано течение газа, зеленым — нефти
Несложно видеть, что соотношения (1) и (2) изначально носят сугубо приближенный характер, так что именно их вид целесообразно устанавливать на основе результатов измерений. Члены и в уравнении (1) соответствуют нестационарности течения и воздействию неоднородного поля давления на флюид соответственно. Очевидно, что они должны входить в уравнение динамики именно в том виде, в каком приведены в уравнении (1). Остальные три члена могут быть сформулированы иначе: хоть и имеет прозрачный физический смысл, но, будучи квадратичным относительно скорости, часто оказывается избыточным из-за своей малости; зависит от угла θ, достоверность значения которого ограничена методологическими погрешностями; является приближенным выражением для силы трения флюида о стенку трубы. Последнее обстоятельство особенно важно: неопределенность в описании трения может привести к появлению в уравнении (1) дополнительных членов произвольной формы, зависящих не только от и , но также и от . Таким образом, в общем случае закон изменения импульса флюида должен выражаться не уравнением (1), а соотношением следующего вида:
. (3)
Здесь F представляет собой неизвестную функцию, к нахождению вида которой и сводится рассматриваемая в данной работе задача. Поскольку желательно, чтобы соотношение (3) в результате реконструкции имело прозрачный физический смысл, входящую в него неизвестную функцию F целесообразно искать в виде линейной комбинации некоторого числа n заранее выбранных функций F1, F2, …, Fn, формирующих тем самым библиотеку допустимых членов:
. (4)
Здесь c1, c2, …, cn представляют собой коэффициенты линейной комбинации, однозначно определяющие вид функции F. Если зафиксировать значения этих коэффициентов, то приведенную выше систему уравнений можно решить численно и далее проверить, насколько близки предсказанные показания имеющихся в скважине датчиков к измеренным величинам (отметим, что на практике точное совпадение не может быть достигнуто из-за конечности размера n библиотеки). Отсюда следует, что изучаемая задача реконструкции может быть сформулирована в терминах минимизации некоторой невязки R, представляющей собой функцию n переменных c1, c2, …, cn. Очевидно, такая численная оптимизация требует расчета значений функции невязки в довольно большом числе различных точек c1, c2, …, cn, для каждой из которых необходимо полностью решить систему определяющих уравнений. Следовательно, для эффективного применения развиваемого подхода как численный решатель системы [6], так и алгоритм оптимизации должны обеспечивать максимально быструю сходимость. Измерения, на основе которых может быть осуществлена идентификация модели, должны относиться к значениям независимых переменных vg, vl, αg, αl и p. При этом первые четыре из них регистрируются расходомером, снабженным датчиками состава, а последнее — датчиками давления. Такой сценарий будем обозначать как «максимум информации». Противоположной является ситуация, когда в скважине имеется единственный датчик давления на забое — этот сценарий условно назовем «минимум информации». Конкретный вид функции невязки должен выбираться в зависимости от набора датчиков, присутствующих в скважине. Так в случае сценария «максимум информации» можно использовать следующее выражение:
. (5)
Здесь и далее знаком «~» обозначены значения, полученные в результате измерений, а суммирование по t ведется по всем моментам времени, в которые эти измерения проводились. Точкам x = xj соответствуют положения m датчиков давления в скважине. Несложно видеть, что невязка (5) представляет собой сумму относительных отклонений (строго говоря, модифицированных, поскольку абсолютное отклонение нормируется не на одну из сравниваемых величин, а на их среднее арифметическое, что позволяет максимально избежать появления знаменателей, близких к нулю) предсказанных показаний. В сценарии «минимум информации» невязка формулируется аналогично.
Рассматриваемый подход является развитием подходов к идентификации моделей, представленных в [7, 8], где предполагалось, что решения уравнений, вид которых подлежит восстановлению, могут быть представлены в виде известных функций как от времени, так и от координаты. Напротив, здесь будет рассмотрен в том числе случай, когда в скважине расположен всего один датчик давления, так что пространственное распределение параметров течения оказывается полностью недоступным для измерений. Возможность идентификации динамической модели в условиях дефицита данных, а также акцент на отсутствии необходимости в установлении тех из уравнений системы, в которых нет неопределенности, составляют методическую новизну данной работы.
Результаты: восстановление уравнений
Чтобы протестировать изложенный подход к реконструкции уравнений, будем использовать в качестве известных результатов измерений численное решение исходной системы. При этом уравнение (1) выполняется точно, так что в случае реконструкции должно получиться
. (6)
Будем использовать следующие значения параметров: L = 1 м, p0 = 105 Па, μ = 0.029 кг/моль, T = 300 К, ρl = 850 кг/м3, g = 9.8 м/с2, θ = 0, k = 10 м-1. Также будем считать, что Qg(t) = 2 кг/(м2∙с), то есть массовый расход газа является постоянным. Относительно массового расхода Ql(t) нефти, напротив, будем предполагать, что он испытывает кратковременное возмущение, так что изначальный расход 1000 кг/(м2∙с) в течение 0.1 с плавно повышается до 1250 кг/(м2∙с) или 1500 кг/(м2∙с), после чего так же плавно за 0.1 с возвращается к исходному значению и далее не меняется. Такой сценарий соответствует управлению добычей месторождения при помощи забойного штуцера. С математической точки зрения наличие двух различных значений максимального расхода связано с тем, что отклик системы на одиночное возмущение не всегда позволяет с достаточной точностью восстановить коэффициенты в формуле (4). Тогда может быть полезно рассмотреть отклик системы на два (или более при необходимости) последовательных возмущения различной амплитуды. При этом фактически сначала осуществляется минимизация невязки для случая одиночного возмущения меньшей амплитуды, причем в качестве начальных значений искомых коэффициентов ci выбираются некоторые произвольные величины (например, ci = 5). Затем полученные оптимальные значения ci уже сами используются в качестве начальных при повторной оптимизации для случая возмущения большей амплитуды. В качестве примера используется простейшая библиотека из n = 4 членов,
, (7)
что приводит к тому, что в сценарии, когда в скважине имеется сетка датчиков давления и расходомер, установленный на устье, за 300 итераций оптимизатора значения коэффициентов c1 = 9.8, c2 = 1, c3 = 10, c4 = 0 восстанавливаются точно даже без привлечения второго возмущения массового расхода нефти. Соответствующий график постепенного изменения коэффициентов ci в процессе оптимизации в зависимости от номера итерации приведен на рисунке 2. В сценарии, когда данные поступают с единственного датчика давления на забое скважины, одного возмущения становится уже недостаточно: за 460 итераций достигаются значения коэффициентов, содержащие заметные погрешности (в частности, c4 = 0.73). Следует, однако, подчеркнуть, что даже при этих значениях предсказанные показания единственного датчика давления оказываются в соответствии с истинными данными. Добавление повторной оптимизации с использованием второго возмущения позволяет за 250 дополнительных итераций полностью избавиться от имеющихся погрешностей, получив верные значения коэффициентов (рис. 3).
 |
 |
|
Рисунок 2. Оптимизация коэффициентов для библиотеки из четырех членов. Сценарий «максимум информации» |
Рисунок 3. Оптимизация коэффициентов для библиотеки из четырех членов. Сценарий «минимум информации» |
Далее рассмотрим расширенную библиотеку из n = 6 членов:
. (8)
В этом случае в сценарии «минимум информации» приемлемый результат не удается получить даже при использовании двух возмущений массового расхода нефти. Тем не менее, в сценарии «максимум информации» оказывается достаточно уже одного возмущения; при этом требуется 800 итераций оптимизатора, что значительно больше, чем было необходимо при выборе библиотеки (7). Графики изменения коэффициентов, аналогичные рисункам 2 и 3, показаны на рисунках 4 и 5. Несложно видеть, что в сценарии «максимум информации» уже примерно за 400 итераций форма уравнения (3) восстанавливается на качественном уровне верно. Следует отметить, что определенный выбор функций-кандидатов для библиотеки может приводить к нефизичным формам законам сохранения, что следует избегать при помощи экспертной оценки.
 |
 |
|
Рисунок 4. Оптимизация коэффициентов для библиотеки из шести членов. Сценарий «максимум информации» |
Рисунок 5. Оптимизация коэффициентов для библиотеки из шести членов. Сценарий «минимум информации» |
Результаты: уточнение параметров
Из смысла рассмотренной задачи идентификации динамической модели ясно, что для фиксированной скважины эта процедура должна выполняться однократно. Потребность в радикальном пересмотре модели может возникнуть лишь в том случае, если условия течения претерпели изменения на качественном уровне, например в результате изменения конструкции скважины. Тем не менее числовые параметры единожды установленной модели могут, вообще говоря, потребовать уточнения даже вследствие чисто количественных изменений в режиме эксплуатации. В первую очередь это связано с тем, что значения ряда параметров, таких как k, C, vd, в действительности не являются строго постоянными и находятся в некоторой зависимости как от скоростей, так и от состава флюида [5]. Следовательно, даже постепенное снижение пластового давления может рано или поздно привести к существенной ошибке в используемых значениях параметров модели, что сделает ее применение некорректным, несмотря на верную форму составляющих модель уравнений. Тем самым оказывается крайне важной задача определения и дальнейшего уточнения значений параметров фиксированной модели в отрыве от задачи ее изначальной идентификации. Очевидно, что данная проблема становится еще более актуальной в тех случаях, когда уравнения модели не восстанавливаются на основе имеющихся измерений, как предлагается в настоящей работе, а постулируются из физических соображений.
Альтернативным подходом, избавленным от указанных недостатков, является частный случай методики, предложенной выше для идентификации динамических систем, и заключающийся в минимизации невязки как функции уточняемых параметров. Чтобы убедиться в этом, будем считать, что вид функции F в формуле (3) уже установлен и совпадает с (1), так что динамика флюида в модели потока дрейфа описывается уравнением (2). Пусть также неизвестными параметрами системы являются k, C, vd. Тогда невязка оказывается функцией этих трех параметров, искомые значения которых могут быть найдены путем численной оптимизации.
 |
 |
|
Рисунок 6. График изменения значений параметров в процессе оптимизации в зависимости от номера итерации для гомогенного течения |
Рисунок 7. График изменения значений параметров в процессе оптимизации в зависимости от номера итерации при наличии проскальзывания фаз |
На рисунках 6, 7 показаны графики изменения значений искомых параметров в процессе оптимизации невязки R для случаев гомогенного течения и наличия проскальзывания между фазами. Результаты приведены для сценария «минимум информации». В обоих случаях были получены верные значения параметров: k = 10 м-1, C = 1, vd = 0 для случая гомогенного течения (рис. 6) и k = 10 м-1, C = 1.2, vd = 0.4 м/с при наличии проскальзывания (рис. 7). Таким образом, задача уточнения параметров модели путем минимизации невязки R решается быстрее, чем задача идентификации модели. Это объясняется как уменьшением числа независимых переменных, по которым осуществляется оптимизация R, так и меньшей изменчивостью уравнений, описывающих течение: ясно, например, что варьирование коэффициента перед членом в уравнении (4) может качественно повлиять на характер получаемых решений, из-за чего вид зависимости R от значений коэффициентов ci существенно усложняется, а численная оптимизация, в свою очередь, становится менее надежной.
Выводы
Предложенный в данной работе подход к идентификации модели физического процесса обладает значительной универсальностью. По сути, он может быть успешно применен к произвольной задаче, для которой физические соображения позволяют хотя бы в общих чертах сконструировать систему, состоящую из уравнений в частных производных и алгебраических уравнений. Далее те из этих уравнений, вид которых неизвестен или известен не полностью, могут быть реконструированы на основе результатов имеющихся измерений. Данную технику можно рассматривать в первую очередь как инструмент физического исследования, позволяющий автоматизировать построение модели, которая заведомо подходит для описания наблюдаемого явления, обладая при этом максимальной простотой и логической обоснованностью. В качестве направлений будущей работы предлагается рассмотреть чувствительность алгоритма к шуму измерений. Кроме этого, необходимо ввести механизм регуляризации, что позволит масштабировать метод на библиотеки больших размеров.
