Стрекалов А.В. Математические модели гидравлических систем для управления системами поддержания пластового давления
Подождите немного. Документ загружается.
251
рых могут принимать нулевые значения в определенных точках (в ча-
стности, в экстремумах). Например, для функции
f(q)=q
3
первая про-
изводная
2
3q)q(f =
′
и при q=0 тоже равна нулю. Естественно СЛАУ
(4.47) в этом случае будет несовместной. Поэтому здесь рекомендуется
конечное выражение для приближенного вычисления производных
i
I
iii
I
ii
i
I
ii
h
qfhqf
q
qf )()()(
)()()(
−+
≈
∂
∂
,
где
h
i
– величина смещения по оси расходов звена i.
Для метода Ньютона часто рекомендуется вычислять величины
h
i
на каждой итерации и для каждой функции по следующей формуле
α
⋅=
ii
qh
.
Согласно опыту автора, такой подход к корректировке смещения
h
i
не только не дает эффекта в большинстве случаев, но и «заводит в ту-
пик» (условия решения не достигаются) процесс сходимости для
большинства моделей гидросистем.
Предлагается «алгоритм динамичного подбора смещения»
h
i
, суть
которого сводится к следующим основным положениям.
1. Смещение
h
i
выбирается одинаковым для всех замыкающих от-
ношений (
f(q) – моделей элементов гидросистемы).
2. При первой попытке свести итерационный процесс к решению
величина
h принимается равной 0.01 м
3
/с (т.е. максимально возможной
для реальных гидросистем). Для массового расхода можно взять 10
кг/с (для жидкостей). Согласно опыту достаточно (для повышения
скорости) взять
h=1⋅10
–6
, чтобы алгоритм долго не подбирал смещение.
3. При достижении критического числа итерации
I > I
max
или при
повышении погрешности на очередной итерации относительно преды-
дущей
)I()I(
εε
>
+1
, h делится на 10 (установлено экспериментально,
что это является лучшим делителем), далее восстанавливается лучшее
из найденных приближений
)(
L
I
i
q для всех звеньев гидросистемы и
процесс решения запускается заново.
4. Для восстановления лучшего приближения ведется история про-
цесса сходимости, т.е. запоминаются на каждой итерации текущие
приближения и соответствующие им относительные (4.50), абсолют-
ные (4.49, 4.51) и интегрированные –
)I(
ε
погрешности, затем после
сортировки массивов по условию увеличения погрешности выбирается
самый первый массив текущих приближений
)0(
i
q и соответствующая
(4.52)
252
ему погрешность. Далее итерационный процесс продолжается с новым
смещением
h от итерации I=1 с восстановленным приближением
)0(
i
q
для всех звеньев.
5. Для учета текущей погрешности величины
r
ϕ
и
j
γ
суммируют-
ся, а затем находятся их средние арифметические:
∑
=
=
c
r
r
1
ϕϕ
,
∑
=
=
t
j
j
1
γγ
.
Причем, по опыту установлено, что погрешность
γ
всегда на несколь-
ко порядков меньше
ϕ
, и составляет для чисел «с плавающей точкой
двойной точности» около
13
101
−
⋅
, что практически соответствует ну-
лю в принятом масштабе используемых величин. Это объясняется тем,
что первая часть системы (4.4) линейна, а погрешность связана лишь с
неточностью метода решения СЛАУ и вычислениями мат. операций с
числами математическим сопроцессором ЭВМ.
На рис. 4.16 показана общая схема описанного алгоритма числен-
ного решения для метода «путевой увязки».
Подытоживая да
нный пункт, следует отметить, что такой подбор
смещения сделал практически универсальным адаптированный метод
Ньютона для решения системы (4.4): на данный момент сбоев в про-
цессе сходимости не обнаружено как для моделей реальных, так и для
моделей сложных экспериментальных гидросистем в программном
комплексе Hydra’Sym[140].
Естественно, описанный выше метод численного решения неодно-
родной системы не
линейных алгебраических уравнений при решения
потокораспределения гидросистем не может работать без предвари-
тельно подготовленной матрицы
B линейно-независимых путей.
253
Рис. 4.16. Схема алгоритма численного решения потокораспределения методом «путевой увязки» (АПУН)
254
Рассмотрим метод численного решения потокораспределения для
условий формирования системы уравнений (4.26) метода «узловой
увязки».
Для решения системы (4.26) можно использовать метод Ньютона,
но как показал опыт многих ученых, данный метод не является уни-
версальным, хотя и при «удачном стечении обстоятельств» показывает
хорошую скорость сходимости и точность.
Методом Ньютона итерационный процесс строится посредством
ли
неаризации функций
F
j
из (4.20) на каждом шаге для текущего при-
ближения давлений транзитивных узлов. На каждом шаге строится
матрица Якобиан (частных производных), которая в системе СЛАУ
умножается на вектор текущих приращений
)(I
j
pΔ
и приравнивается к
вектору текущих значений функций
F
j
. Ниже приводится СЛАУ отно-
сительно приращения
)I(
j
p
1+
Δ
для итерации I+1 для решения системы
(4.20).
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎧
=Δ
∂
∂
+
++Δ
∂
∂
++Δ
∂
∂
=Δ
∂
∂
+
++Δ
∂
∂
++Δ
∂
∂
=Δ
∂
∂
+
++Δ
∂
∂
++Δ
∂
∂
+
++
+
++
+
++
)p,p,..,p,p(Fp
p
)p,p,..,p,p(F
..p
p
)p,p,..,p,p(F
..p
p
)p,p,..,p,p(F
)p,p,..,p,p(Fp
p
)p,p,..,p,p(F
..p
p
)p,p,..,p,p(F
..p
p
)p,p,..,p,p(F
)p,p,..,p,p(Fp
p
)p,p,..,p,p(F
..p
p
)p,p,..,p,p(F
..p
p
)p,p,..,p,p(F
)I(
m
)I(
j
)I()I(
t
)I(
t
t
)I(
m
)I(
j
)I()I(
t
)I(
j
j
)I(
m
)I(
j
)I()I(
t
)I(
)I(
m
)I(
j
)I()I(
t
)I(
m
)I(
j
)I()I(
j
)I(
t
t
)I(
m
)I(
j
)I()I(
j
)I(
j
j
)I(
m
)I(
j
)I()I(
j
)I(
)I(
m
)I(
j
)I()I(
j
)I(
m
)I(
j
)I()I()I(
t
t
)I(
m
)I(
j
)I()I(
)I(
j
j
)I(
m
)I(
j
)I()I(
)I(
)I(
m
)I(
j
)I()I(
21
1
21
1
21
1
1
1
21
21
1
21
1
21
1
1
1
21
211
1
211
1
211
1
1
1
211
M
M
Таким образом, для последовательных решений (4.52) начальные
значения давлений
)0(
j
p к первому – I = 0 приближению должны быть
(4.52)
255
заданы, а последующие
)I(
j
p
будут рассчитываться. Величины давле-
ний в активных узлах являются константами на всех итерациях.
Характер способа нахождения частных производных фактически
полностью определяет скорость и качество сходимости метода Ньюто-
на. Здесь частные производные лучше находить в конечно-разностном
виде
()
I
j
)I(
m
)I(
j
)I()I(
j
)I(
m
)I(
j
)I(
j
)I()I(
j
j
)I(
m
)I(
j
)I()I(
j
h
)p,p,..,p,p(F)p,hp,..,p,p(F
p
)p,p,..,p,p(F
2121
21
−+
≈
≈
∂
∂
,
однако, в отличие от метода «путевой увязки», это далеко не всегда
приводит к инвариантной сходимости.
Решив (4.52) любым методом решения СЛАУ, получаем прираще-
ния к текущим приближениям неизвестных, а затем их складываем:
)()()1( I
j
I
j
I
j
ppp Δ+=
+
.
После переходим к формированию новой СЛАУ – (
I+1).
Так до тех пор, пока не выполнятся условия точности: абсолютного
материального баланса в транзитивных узлах
tjabsQpS
ji
ii
,...2,1,)( =<Δ
∑
∈
,
и относительного материального баланса
t,...,j,
)p(S
)p(S
ji
ii
ji
ii
21=<
Δ
Δ
∑
∑
∈
∈
ε
где
absQ – величина расхода, принимаемая в качестве допущения,
как ноль (абсолютная погрешность);
ε
– допустимая относи-
тельная погрешность;
i
p
Δ
– перепады давлений в звеньях, смеж-
ных с узлом
j.
Здесь также следует отметить, что выбор величины
)(I
j
h
для вычис-
ления частных производных определяет качество сходимости. Вычис-
ление величины
)(I
j
h
на каждой итерации по следующей формуле
ε
⋅=
)()( I
j
I
j
ph
,
(4.53)
(4.54)
(4.55)
(4.56)
(4.57)
256
как и для предыдущего метода в большинстве случаев эффекта не да-
ет.
По аналогии с рассмотренным ранее методом здесь предлагается
динамичный подбор смещения
h,
суть которого сводится к следующим
основным положениям.
1. Смещение
h выбирается одинаковым для всех функций F
j
.
2. Вначале
h принимается равным 10000 Па.
3. При достижении критического числа итерации
I>I
max
или при по-
вышении погрешности на очередной итерации относительно преды-
дущей
h делится на 10, далее восстанавливается лучшее из найденных
приближений
)(
L
I
j
p
для всех неизвестных давлений.
4. Для восстановления лучшего приближения ведется история про-
цесса сходимости, т.е. запоминаются на каждой итерации текущие
приближения и соответствующие им относительные погрешности, по
достижении условия 3 итерационный процесс продолжается с новым
смещением
h
(I+1)
=h
(I)
/10 и с восстановленным лучшим приближением
)0(
j
p
.
5. Для учета текущей погрешности вводятся величины текущих от-
носительных погрешностей, найденного баланса потоков в узлах
t,...,j,
)p(S
)p(S
e
ji
ii
ji
ii
j
21=
Δ
Δ
=
∑
∑
∈
∈
,
которые суммируются, а затем находится их среднее арифметическое
t
e
e
t
j
j
∑
=
=
1
, показывающее интегрированную относительную погреш-
ность текущей итерации.
При всей кажущейся простоте описанного метода решения далеко
не всегда можно получить хорошо сходящийся или вообще сходящий-
ся к решению итерационный процесс.
Далее рассмотрим наиболее эффективные варианты адаптации
данного метода к решению СНАУ (4.26).
Важным моментом для расчета потокораспределения является опи-
сание замыкающих отношений, входящих в ту ил
и иную формы запи-
си системы уравнений потокораспределения.
(4.58)
257
Теоретический расчет (моделирование) любого объекта (элемента)
гидросистемы сводится к определению перепада давления между вхо-
дом в объект и выходом из него в зависимости от внутренних свойств
объекта. Иначе говоря, целью такого расчета является нахождение за-
висимости между величиной расхода текучей среды сквозь объект
q
i
и
перепадом давления на его концах, зависящим от свойств объекта –
)(
iii
qfp =
′
Δ
. Перепад давления
i
p
′
Δ
, зависящий от внутренних
свойств объектов, это та часть общего перепада давления, которая име-
ет зависимость от расхода. Остальная составляющая
Δz
i
перепада Δp
i
не имеет зависимости от расхода, поэтому далее перепад
i
p
′
Δ
будем
обозначать просто
i
pΔ
. Итак,
i
pΔ
есть разность между давлением на
входе в звено и давлением на выходе из него, т. е. знак
i
pΔ
строго за-
висит от ориентации объекта в структуре гидросистемы, его влияния
на энергию потока текучей среды и направления течения потока.
К виду функций
)(
ii
qf во всех существующих методах численно-
го решения потокораспределения предъявляются достаточно жесткие
требования.
1.
0)0( =
i
f
.
2. Первая производная:
0≠
′
)q(f
i
на всей области определения
функции.
3.
)()( qfqf
ii
−=−
, т.е. должна быть нечетной.
Такие условия невозможно соблюсти для большинства моделируе-
мых элементов гидросистем: труб (в автомодельной области течения,
насосов на всех режимах, скважин, запорной арматуры, штуцеров, зон
пласта с квадратичным законом фильтрации и т.п.).
В связи с этим в большинстве существующих моделей это ограни-
чение обходится фиктивным учетом та
ких объектов: описание насосов
фиксированным давлением, скважин – постоянным оттоком или при-
током текучей среды. Все эти допущения уводят нас от главной цели
адекватной имитации гидросистемы.
Такие принимаемые учеными допущения связаны, прежде всего, с
недостаточной универсальностью математических методов решения
систем нелинейных уравнений.
Как уже отмечалось, основным методом решения систем нелиней-
ных ал
гебраических уравнений является метод Ньютона, который
имеет ряд существенных недостатков.
258
1. Сходимость метода определяется точностью начального при-
ближения, т.е. при недостаточно точном (трудно заранее сказать, ка-
кой точности будет достаточно) определении начального приближения
метод даже «не пытается» сходиться.
2. Линеаризация функций, имеющих экстремумы в районе текуще-
го приближения, приводит к бесконечному итерационному процессу.
3. Для надежной сходимости метода необходим достаточно то
чный
алгоритм выбора шага
h при численном определении частных произ-
водных замыкающих отношений в методе путевой увязки
h
qfhqf
qf
iiii
ii
)()(
)(
−+
=
′
и функций
=
′
),,..,,(
)()()(
2
)(
1
I
m
I
j
II
j
ppppF
=
h
ppppFphpppF
I
m
I
j
II
j
I
m
I
j
II
j
),,..,,(),,..,,(
)()()(
2
)(
1
)()()(
2
)(
1
−+
в методе узловой
увязки.
4. Невозможно получить сходящийся к решению процесс, если су-
ществует большой разброс условных проводимостей звеньев.
5. Невозможно получить сходящийся к решению процесс, если в
структуре гидросистемы имеются элементы с замыкающими отноше-
ниями, функции которых по-разному ведут себя при
q<0 и q>0, напри-
мер, обратные клапаны.
6. Несовместимо с решением наличие закрытых задвижек (со
сверхбольшим сопротивлением) на выкидной линии насосных агрега-
тов.
Рис. 4.17. Типовые гидравлические характеристики (замыкающие отно-
шения) элементов гидросистем: а – трубопровод; б – элемент пласта
(ячейка);
в – обратный клапан
259
Рассмотрим наиболее распространенные формы гидравлических
характеристик (замыкающих отношений) элементов гидросистем (рис.
4.17). Как видно из рис. 4.17, функции гидравлических характеристик
по большей части нелинейны, а например, для трубопроводов –
0)0( =
тр
f
. Особенно сложными для решения являются характеристи-
ки обратных клапанов (рис. 4.17 –
в), знак первых производных кото-
рых изменяется многократно в районе
q = 0.
Как видно из системы (4.20), при ее решении необходим поиск
функций
(
)
tjj
ppppF ,..,,..,,
21
, сложенных функциями
(
)
ii
pS
Δ
или
)z,z,p,p(S
ieibieib
jjjji
′
, которые, в свою очередь, являются решением
уравнений
()
0
=
Δ−
iii
pqf
относительно q
i
в том или ином виде.
Для выявления характера функций
(
)
tjj
ppppF ,..,,..,,
21
построим (из
примера на рис. 4.7/2) кривую функции
()
65422
,,, ppppF
при измене-
нии
p
2
и постоянных остальных давлениях: p
5
= const =0, p
4
= const =0,
p
6
= const =0 (см. рис. 4.18).
Данный распространенный вид для функций
F
j
(p…) вызывает у хо-
рошо известного метода Ньютона трудности в сходимости при произ-
вольном начальном приближении величин давлений транзитивных
Рис. 4.18. Гидравлическая характеристика узла 2
260
узлов. На рис. 4.18 также показан пример проблемы сходимости мето-
да Ньютона для подобных функций при нулевом начальном прибли-
жении, из которого видно, что при попадании текущего приближения
в зоны кривой с близкими или равными значениями первой производ-
ной решение не может быть достигнуто. В дополнение к вышеописан-
ному недостатку метода Ньюто
на следует также отметить, что при
нахождении частных производных локально, т.е. в соответствии с ве-
личиной конечного приращения аргумента
h или алгебраически (что
не всегда возможно), текущее приближение неизвестного к корню,
получаемое на их основе, может оказаться дальше, чем предыдущее.
В связи с вышесказанным для решения системы (4.26) требуются
выбор и адаптация известных методов решения СНАУ к условиям дан-
ной задачи потокораспределения.
Ниже предлагаются два метода решения систем нелинейных алгеб-
раических уравнений, к
оторые позволят снять большинство ограниче-
ний на вид замыкающих отношений
)(
ii
qf
и структурных особенно-
стей гидросистем.
Здесь предлагаются модификации (или адаптации) метода Ньюто-
на, названные фрагментарной линеаризацией, незначительно повы-
шающей надежность сходимости метода Ньютона, и адаптивной ли-
неаризацией, являющиеся практически панацеей для решения систем
нелинейных уравнений в плане надежности и скорости сходимости
итерационного процесса последовательных приближений.
Рассмотрим способ фрагментарной линеаризации, суть которого
своди
тся к предварительной обработке исходных замыкающих отно-
шений
()
qf
или
(
)
pS Δ
(в зависимости от исходной системы уравне-
ний) таким образом, чтобы получить достаточно четкую сходимость.
Для этого итоговые функции должны:
1. быть линеаризованы до начала итерационного процесса;
2. в областях определения аргумента, где заведомо не может быть
корня, первые производные (касательные) должны быть направлены в
начало координат;
3. в зонах четности первой производной ис
ходной гидравлической
характеристики
f(q) –
)()(
RiLi
qqfqqf >
′
=<
′
или производной об-
ратной ей функции
S(Δp) –
(
)
(
)
RL
ppSppS
Δ
>
Δ
′
=
Δ
<
Δ
′
линеаризо-
ванные функции должны иметь различные величины производных:
искусственно создаются флуктуации производных во избежание про-