Стрекалов А.В. Математические модели гидравлических систем для управления системами поддержания пластового давления
Подождите немного. Документ загружается.
261
блемы на рис. 4.18, где четность производной проявляется при
L
p = –
5,
R
p =15.
На рис. 4.19 показан пример фрагментарной линеаризации функции
)(
33
pS Δ
.
Как видно из рис. 4.19, вне зоны приемлемых корней касательные
направлены в начало координат, а в зонах четности первых производ-
ных –
()
(
)
RL
ppSppS
Δ
>
Δ
′
=
Δ<Δ
′
касательные имеют разные углы.
В зависимости от требуемой точности функции замыкающих от-
ношений могут быть линеаризованы более или менее плотными фраг-
ментами. Здесь детально не рассматриваются вопросы практического
применения данного способа в связи с тем, что алгоритмизация на-
столько сложна (это свя-
зано с наличием множе-
ства логических проце-
дур), что рассмат
ривать
ее в данной работе не
имеет смысла.
Безусловно, рассмот-
ренный способ снижает
точность решения, од-
нако более стабилен в
плане сходимости.
Более удачным спо-
собом адаптации являет-
ся адаптивная линеари-
зация. Суть метода сво-
дится к использованию
при расчете в методе
«путевой увязки» вместо частных производных
i
)I(
ii
q
)q(f
∂
∂
их конечные
выражения вида
)I(
)I()I(
ii
)I()I(
ii
H
)Hq(f)Hq(f
2
−−+
,
где
H
(I)
– величина, задаваемая для всех неизвестных константой
(для текущей итерации), которая охватывает область определе-
ния всех неизвестных в районе текущего приближения каждой.
Фрагмент
Рис. 4.19. Пример фрагментарной линеаризации
262
Вместо функций
)(
ii
qf
на каждой итерации I как при численном
решении СНАУ метода «путевой увязки», так и при решении (4.26)
методом «уловой увязки» используется их линейное преобразование
)(
)(I
i
L
i
qf
:
)(fq
H
)Hq(f)Hq(f
)q(f
i
)I(
i
)I(
)I()I(
ii
)I()I(
ii
)I(
i
L
i
0
2
+
−−+
=
или для условия течения сжимаемой ТС
)(fM
H
)HM(f)HM(f
)M(f
i
)I(
i
)I(
)I()I(
ii
)I()I(
ii
)I(
i
L
i
0
2
+
−−+
=
.
Причем
DHH
II
/
)()1(
=
+
, где D – параметр сходимости, который
чем больше, тем меньше скорость сходимости, но больше стабиль-
ность. Оптимальными значениями
D с точки зрения автора является
]3,5.1[∈D
.
Если в итерационном процессе
I>I
max
и не достигается требуемая
точность, то восстанавливается лучшее приближение
)0(
i
q или
)0(
j
p и
величина
)(I
H
корректируется с меньшим, чем ранее (в итерации с
наилучшим запомненным приближением) значением
D.
Для решения задачи (4.26) используется тот же прием с функциями
)(qf
i
. Таким образом, они линеаризуются и уже, будучи линеаризо-
ванными, используются на следующих этапах вычислений текущей
итерации в виде (4.59).
Применение такого рода линеаризации позволяет избежать посто-
янного решения (на каждой итерации) нелинейных уравнений
(
)
0=Δ− pqf
для нахождения функций
(
)
pSq
Δ
=
, они будут опи-
сываться линейным соотношением
)Hq(f)Hq(f
H)(fp
)p(S
)I()I(
ii
)I()I(
ii
)I(
ii
i
)I(
i
−−+
−Δ
=Δ
2
1
0
или для условия течения сжимаемой ТС
)HM(f)HM(f
H)(fp
)p(S
)I()I(
ii
)I()I(
ii
)I(
ii
i
)I(
i
−−+
−Δ
=Δ
2
1
0
.
Такого рода прием относится и к нахождению функций
)z,z,p,p(S
ieibieib
jjjji
′
, вычисление которых основано на расчетах f
i
(q),
и, следовательно, подразумевает замену их на (4.59).
(4.59)
(4.60)
263
Особенно важно, что для решения задачи (4.26) после нахождения
нового приближения необходимо рассчитать соответствующие ему
значения расходов
q
i
для того, чтобы провести линеаризацию (4.59) и
(4.60) к началу следующей итерации.
Таким образом, для решения (4.26) адаптивная линеаризация необ-
ходима лишь для основных замыкающих отношений, т.е. гидравличе-
ских характеристик
)(qf
i
. Начальное приближение данного способа
метода (4.52) лучше выбрать нулевым, т.е.
0
)0(
=
j
p
.
На рис. 4.20 показан пример решения одного нелинейного уравне-
ния
f(q)=0 с использованием адаптивной линеаризации с параметром
D=3.
Метод адаптивной линеаризации сочетает в себе высокую скорость
сходимости метода Ньютона и высокую надежность сходимости при
нулевом первом приближении аналогично методу деления отрезка по-
полам («бисекций») с исключением требования к охвату области корня
величиной
H
(I)
.
На рис. 4.21 приведена схема алгоритма численного решения пото-
кораспределения для условий течения сжимаемой ТС с применением
Рис. 4.20. Сходимость метода адаптивной линеаризации для при-
мера решения нелинейного уравнения
f (q) =0
264
адаптивной линеаризации.
Во избежание некоторых разночтений следует отметить, что при
расчете потокораспределения в условиях течения сжимаемой жидко-
сти под замыкающими отношениями
f
i
(M
i
) подразумевается
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
i
i
i
M
f
ρ
,
если по звену не изменяется плотность и
∑
=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
Δ
N
k
i
i
i
l
M
f
k
1
,
ρ
или
)(
ii
Mw , полученные из (4.28), если учитывается изменение объемно-
го расхода и плотности по длине звена. Для неизотермического тече-
ния данные замыкающие отношения вычисляются исходя из известно-
го на данном шаге решения распределения температуры потока в каж-
дом звене (см. далее).
При построении СЛАУ в методе узловой увязки алгоритма
«АЛУС» необходимо найти значения част
ных производных –
=
∂
∂
j
I
m
I
j
II
j
p
ppppF ),,..,,(
)()()(
2
)(
1
=
h
)p,p,..,p,p(F)p,hp,..,p,p(F
)I(
m
)I(
j
)I()I(
j
)I(
m
)I(
j
)I()I(
j 2121
−+
,
где значение
h для данного алгоритма особого значения не имеет и
может быть рассчитано по формуле
h = H
(I)
/100.
Для расчета значений самих функций
(
)
tjj
ppppF ,..,,..,,
21
необхо-
димо рассчитать соответствующие им функции
),,,(
ieibieib
jjjji
zzppS
или
)z,z,p,p(S
ieibieib
jjjji
′
[см. выражения (4.20а) и (4.21)] . Для их рас-
чета на каждой итерации происходит линеаризация функций
f
i
(M) та-
ким образом, что звено
i разделяется на N частей и полагается, что
f
i
(q,Δl) описывают перепад давления на каждом участке k на длину
l
Δ
при постоянстве плотности и вязкости [см. формулу (4.28)]. Причем,
на каждой итерации изменяется параметр линеаризации –
H. При дос-
тижении критического числа итераций параметр
D для модификации
H снижается в два раза, и процесс запускается с наилучшего прибли-
жения.
265
Рис. 4.21. Схема алгоритма адаптивной линеаризации для метода «узловой увязки» условия течения
сжимаемой ТС (АЛУС)
266
В общем смысле (если не использовать адаптивную линеаризацию)
для нахождения
)z,z,p,p(S
ieibieib
jjjji
′
по известной функции w
i
(M
i
) не-
обходимо решить нелинейное уравнение
)M(wpp
iijj
ieib
=
−
относи-
тельно
M
i
при заданных к текущему приближению давлениях
ib
j
)I(
p
и
ie
j
)I(
p
. Численное решение таких уравнений можно провести методом
Ньютона для решения нелинейных уравнений с одним неизвестным. В
методе «АЛУС» такое решение сводится к решению линейного урав-
нения в виде
)HM(f)HM(f
H
)(fp
)p(S
)I()I(
ii
)I()I(
ii
)I(
ii
i
)I(
i
−−+
−Δ
=Δ
2
1
0
,
где
ieib
j
)I(
j
)I(
i
ppp −=Δ
, а
)p(S
i
)I(
i
Δ
может быть заменена функцией
)z,z,p,p(S
ieibieib
jjjji
′
для течения сжимаемой ТС.
Модель неизотермического течения
Как уже отмечалось, суть расчета комплексного потокораспределе-
ния неизотермического течения сводится к последовательному расчету
потокораспределения в виде решения (4.4) алгоритмом «АПУН» или
(4.22, 4.24, 4.27) алгоритмом «АЛУС» при заданном текущем прибли-
жении термораспределения. Расчет потокораспределения повторяется
до тех пор, пока не выполнится условие точности: соответствие пото-
кораспределения и термораспределения условиям распределения
свойств ТС во всех зв
еньях модели, зависящих от температуры и дав-
ления –
),( Tp
ρ
и
),( Tp
ν
.
Рассмотрим алгоритм численного решения комплексного потоко-
распределения с условием неизотермического течения сжимаемой ТС.
Изначально заданными величинами являются температуры в ак-
тивных узлах, где в результате расчета потокораспределения к теку-
щему приближению может быть выявлен приток.
Процесс последовательных приближений комплексного потокорас-
пределения строится также итеративно. На каждом шаге
J пересчиты-
вается потокораспределение, на основании которого также методом
последовательных приближений определяется распределение темпера-
туры.
267
На текущем шаге приближений
J распределения температуры в ка-
честве основных неизвестных будем полагать глобальное распределе-
ние температуры в транзитивных и активных (где в результате расчета
потокораспределения к текущему приближению –
J выявлен отток из
гидросистемы) узлах –
)J(
j
T
.
Для расчета термораспределения, соответствующего найденному
потокораспределению, необходимо методом последовательных при-
ближений на итерации K+1 найти распределение температур –
)K)(J(
j
T
1+
относительно текущего к шагу J потокораспределения и
предыдущего термораспределения
)K)(J(
j
T
.
Согласно (4.43) температура в узле j на шаге K+1 будет описывать-
ся соотношением метода простой итерации
)T,...,T,...,T,T(T
)K)(J(
m
)K)(J(
j
)K)(J()K)(J(
j
)K)(J(
j 21
1
τ
=
+
.
Итерационный процесс завершается, если для всех искомых темпе-
ратур выполняется условие
absTTT
)K)(J(
j
)K)(J(
j
<−
+1
,
где absT – абсолютная требуемая точность температуры, напри-
мер, 0.001
о
С.
Таким образом, расчет термораспределения (4.61) на каждом шаге
–
J описывается итерационным алгоритмом, вложенным в основной
алгоритм комплексного потокораспределения. После такого распреде-
ления температуры, соответствующего потокораспределению
J, опре-
деляется распределение температуры по звеньям в соответствии с
(4.33) и со всеми вытекающими зависимостями.
Далее к следующему шагу J+1 имеем известным термораспределе-
ние
)J(
j
T
1+
, согласно которому определяем новое потокораспределе-
ние. Процесс повторяется до тех пор, пока не выполнится условие для
каждого узла
absTTT
)J(
j
)J(
j
<−
+1
и для всех узлов
STTT
m
j
)J(
j
)J(
j
<−
∑
=
+
1
1
,
(4.61)
(4.62)
(4.63)
(4.64)
268
Рис. 4.22. Перемерзание потока в звене
где ST – суммарная абсолютная точность по всей системе, например,
0.05
о
С.
Необходимо также отметить, что после нахождения термораспре-
деления может обнаруживаться факт замерзания ТС в любом из звень-
ев, т.е.
кр
)J(
j
jiji
)J(
j
T)T,M(T
ii
<+
∈∈
00
θ
,
где
)J(
j
i
T
0
– температура на входе в звено;
кр
T
– температура замерзания ТС.
Если условие (4.65) выполняется для звена i, то в этом звене зада-
ется кинематическая вязкость
ν
i
=10000 м
2
/с (что соответствует прак-
тически нетекучей среде),
процесс расчета переходит
потокораспределению. Для
описания перемерзших
звеньев следует различать
температуру в узлах начала
и конца от температуры в
начале звена и его конце.
Это связано с тем, что звено
может быть замерзшим, а в
его узлах могут быть
подключены другие звенья,
в которых пе
ремерзания
потока может и не быть (см.
рис. 4.22).
На рис. 4.23 показана схема алгоритма численного решения задачи
комплексного потокораспределения. Параметр
Ф (начальное прибли-
жение искомых температур) может быть установлен в пределах от
кр
T
до
Т
кипения
, например Ф=10
о
С.
Как видно из схемы данного алгоритма, при обнаружении факта
перемерзания в звене алгоритм проверяет, было ли обнаружено пере-
мерзание на предыдущих приближениях, и если такового не было, то
алгоритм изменяет статус только этого звена, не проверяя остальные.
Это связано с тем, что возможно перемерзание этого звена могло из-
менить картину п
отокораспределения, которая может показать на сле-
дующем приближении совершенно другое термораспределение.
(4.65)
269
Например, если один из трубопроводов в коллекторе перемерзает
первым, то в результате потокораспределения окажется, что скорости
потоков в остальных трубопроводах коллектора выросли, что снизит
вероятность их замерзания.
Для расчета динамической модели с учетом изменения внутренних
свойств отдельных или всех элементов ГС алгоритм численного реше-
ния на примере учета регуляторов расхода и д
авления строится со-
гласно следующей схеме.
Рис. 4.23. Алгоритм «АКПР» комплексного потокораспределения (КПР)
270
Схема алгоритма динамической модели с распределенными
параметрами на примере учета регуляторов расхода
Как видно из схемы алгоритма, расчет состояния регулирующих
устройств УП
σ
сводится к последовательным приближениям КПР с
учетом изменения фактора внутреннего состояния регулятора от рас-
хода жидкости, установившегося к данному приближению –
)q(
ii
ζ
.
Момент завершения алгоритма описывается условием абсолютной
погрешности
ζ
abs
определения фактора
i
ζ
. Значение
i
abs
ζ
может
быть выбрано в соответствии с требуемой точностью и условий рабо-
ты каждого звена–регулятора i. Например, для регулятора расхода,
работающего посредством воздействия на поток ТС дросселировани-
ем, величина
i
abs
ζ
может быть принята в 0.00001 м для фактора
i
ζ
–
положение затвора или диаметр диафрагмы.
Численное решение динамической модели ТГС с учетом взаимо-
связи с моделью ГПП рассматривается в разделе 4.7.