Rybinsk, Russian Federation
Rybinsk, Russian Federation
GRNTI 55.01 Общие вопросы машиностроения
GRNTI 55.13 Технология машиностроения
The aim of the work consists in the reveal of the most stable method of identification, in the reveal of factors influencing the parameter definition result of the melt model during leucosapphire singlecrystal growth by Kyropoulos’ method. Investigation methods are based on an engineering process computer modeling under the assumption that the object of control is an aperiodic link of the second order. The comparison of identification results is carried out through different methods built-in Matlab system, their adequacy is defined by actual processes. As a result of investigations there is obtained a description of an aperiodic link of the second order in the form of states space, an interpretation of results obtained in the form of a time constant of object control is given. It is shown that the engineering process description obtained in the form of states space is adequate, the availability of non-zero initial conditions (non-zero speed of technological parameter changes) has no influence upon an identification process. There are ob-tained dynamic characteristics for a set of experimental data by Process Models and Space State methods. The application of Space State method is recommended.
identification, states space, initial conditions, leucosapphire, convective heat exchange, Matlab
Введение
При выращивании монокристаллов лейкосапфира методом Киропулоса типовой технологический процесс выращивания проходит несколько этапов: получение расплава (разогрев), затравление, рост и охлаждение. Затравление и начало роста являются самыми ответственными этапами. Здесь часто возникают ошибки, связанные либо с очень большим сбросом мощности оператором, ведущим к залипанию кристалла, либо с очень малым, ведущим к остановке роста и даже отрыву кристалла от расплава [1–3]. Конструкция ростовой установки и метод выращивания не дают возможности непосредственного измерения температуры. Поэтому на первый план выступают косвенные методы, например по изменению показаний измерителя мощности. По характеру изменения мощности на нагревателе можно определить закон изменения температуры [4]. То есть появляется возможность по результатам косвенных измерений управлять процессом выращивания, например используя модель в составе беспоисковой адаптивной системы стабилизации скорости кристаллизации [5].
Для автоматизации данного процесса решаются различные задачи, одной из которых является получение динамических характеристик объекта управления. Причем идентификация в общем случае может быть проведена для неизвестной формы описания модели объекта [6; 7]. В нашем случае с учетом тепловой природы процесса примем, что модель расплава представляет собой апериодическое звено второго порядка. Проведем идентификацию постоянных времени объекта с помощью пакета ident системы инженерных расчетов Matlab [8] по данным напряжения и мощности нагревателя. Применение способа идентификации Process Models в некоторых случаях даёт неудовлетворительные результаты, поэтому было приято решение использовать для идентификации аппарат пространства состояний, а также проверить влияние начальных условий (НУ) на результаты идентификации.
Математическое описание пространства состояний
В современных пакетах математической обработки форма пространства состояний базируется на матричном представлении и представляет собой наиболее удобную и компактную форму описания многозонных объектов. В пространстве состояний математическое описание системы задается в векторно-матричной форме [9]:
, (1)
где X – вектор состояния размера n×1; U – вектор входных воздействий размера m×1;
A – матрица системы размера n×n; B – матрица входа размера n×m.
В то же времят данная форма дает лишь универсальный механизм, конкретное содержание которого зависит от природы описываемого объекта. Рассмотрим представление в форме пространства состояний объекта второго порядка с передаточной функцией
, (2)
где t1, t2 – постоянные времени объекта; Kp – коэффициент передачи объекта; Tp1 = t1t2; Tp2 = t1+t2.
Такому объекту соответствует дифференциальное уравнение вида
. (3)
Поскольку кроме текущего состояния объекта необходимо учитывать и скорость его изменения, примем вектор состояния в виде
,
где P – мощность, выделяемая на нагревателе, Вт; – скорость ее изменения, Вт/с.
Тогда уравнение для пространства состояний примет вид
. 4)
Определим размерности матриц, входящих в это выражение. Вектор-столбец в левой части имеет размерность 2×1, т.е. для сохранения размерности матрица A должна иметь размерность 2×2. Поскольку входное воздействие (напряжение на нагревателе) является одномерной величиной, то матрица B должна иметь размерность 2×1. Выразим эти матрицы в общем виде через их элементы с учетом соглашения системы Matlab [10] (первый индекс кодирует номер строки, а второй – номер столбца):
,
.
С учетом принятых обозначений представим (4) в виде
. (5)
Матричное уравнение (5) можно преобразовать к системе дифференциальных уравнений в соответствии с правилами перемножения матриц:
(6)
Определим значения коэффициентов матриц A и B таким образом, чтобы получить структуру уравнения (3). Поскольку уравнение (3) второго порядка, то в системе (6) нужно оставить только второе уравнение. Для этого первое уравнение в ней сведем к тождеству, т.е. положим a11 = 0, a12 = 1, b1 = 0. Уравнение (2) приведем к виду
,
а также перенесем в правую часть все, кроме выражения для второй производной
.
Тогда возможно определить недостающие коэффициенты матриц A и B. В итоге они примут вид
,
,
, D = 0.
Идентификация с помощью пространства состояний, влияние НУ
Прежде всего проверим адекватность полученной модели. Для этого сравним реакцию модели объекта второго порядка (2), созданную с помощью функции tf, с реакцией модели в форме пространства состояний (State Space). На вход того и другого объекта подадим линейно нарастающий сигнал с помощью функции lsim. Зададимся постоянными времени t1 = 2500 с, t2 = 50 с, коэффициентом передачи Kp = 5. Для случая нулевых НУ графики переходных процессов совпадают (рис. 1).
Если же задать ненулевую начальную скорость, например 0,01 Вт/с, то начальная стадия переходного процесса существенно изменится (рис. 1). Также меняется вид графика скорости изменения мощности (рис. 2). В системе Matlab начальные условия задаются с помощью функции initial.
После экспорта в рабочее пространство Workspace и преобразования вида tf(n4s2) будет получена передаточная функция
.
После приведения ее к каноническому виду (1 при S в нулевой степени в знаменателе) получим, что
Kp = 5, Tp1 = 124790, Tp2 = 2548.
Учитывая, что t1 = Tp1/ Tp2, t2 = Tp2 - t1, можно утверждать, что
t1 = 48,978 c, t2 = 2499 c.
Аналогичным образом проведем идентификацию с ненулевыми и нулевыми НУ. Результаты приведены в табл. 1, где FTE – процент соответствия полученной модели исходным данным.
Таблица 1
Результаты идентификации
|
НУ |
Kp |
t1, с |
t2, с |
FTE, % |
|
Нулевые |
5 |
48,978 |
2499 |
100 |
|
Ненулевые |
5 |
48,978 |
2499 |
100 |
По данным табл. 1 можно сделать вывод, что наличие НУ не влияет на результаты идентификации параметров объекта. Таким образом, при идентификации реального объекта НУ можно не учитывать.
Идентификация экспериментальных данных
На рис. 4 показан фрагмент экспериментальных данных для напряжения и мощности нагревателя, полученных при выращивании монокристалла. Характерной особенностью данных зависимостей является то, что мощность меняет своё значение с 42,5 до 42,8 кВт, в то время как напряжение поддерживается на уровне 8,15 В (этому предшествовал линейный сброс напряжения). Данный эффект можно объяснить большой инерционностью конвективных потоков в расплаве, которые выносят тепло со дна тигля на поверхность расплава, изменяя распределение температур по высоте вольфрамового высокотемпературного нагревателя. Процесс выращивания монокристалла чувствителен к точности поддержания температуры, и отклонение её значения от поддерживаемой на величину более 0,1 °C является критичным. Поскольку из-за высокой температуры расплава её измерение напрямую невозможно, то используется косвенный метод на основе данных мощности, изменение которой примерно на 60 Вт соответствует 1°C.
Кроме приведенного на рис. 4, рассматриваемый процесс выращивания монокристалла содержит 5 характерных участков со сходной зависимостью мощности от напряжения. Для оценки стабильности методов идентификации Process Models и Space State проведем обработку всех пяти участков, а полученные данные сведём в табл. 2 и 3 соответственно. При анализе результатов будем учитывать не только процент соответствия полученной модели исходным данным (FTE), но и разброс получаемых значений, а также ожидаемый исходя из теплофизических характеристик ростовой установки диапазон постоянных времени (для t1 – десятки секунд, для t2 – тысячи секунд).
Таблица 2
Результаты идентификации с помощью Process Models
|
Участок |
Kp |
|
|
FTE, % |
|
1 |
2,12 |
277,190 |
22906 |
96,66 |
|
2 |
3,22 |
0,007 |
7011 |
98,48 |
|
3 |
5,22 |
0,918 |
10-6 |
36,11 |
|
4 |
5,11 |
0,028 |
1640 |
74,28 |
|
5 |
4,93 |
4,697 |
1169 |
41,42 |
Таблица 3
Результаты идентификации с помощью State Space
|
Участок |
Kp |
|
|
FTE, % |
|
1 |
2,14 |
6,772 |
23445 |
98,02 |
|
2 |
3,23 |
37,408 |
6781 |
97,75 |
|
3 |
5,23 |
52,574 |
1490 |
54,49 |
|
4 |
5,12 |
77,175 |
1362 |
71,60 |
|
5 |
4,93 |
78,647 |
1438 |
40,23 |
Заключение
Сравнение значений t2 в табл. 2 и 3 показывает, что идентификация при помощи Space State позволяет получить более достоверные результаты, о чём свидетельствует отсутствие значений, резко отличающихся от значений на соседних участках. Также получены более корректные значения t1, поскольку предполагается, что она равна десяткам секунд. Полученная модель позволяет косвенно оценить интенсивность движения тепловых потоков, создающих конвективный теплообмен, что может быть использовано в качестве корректирующего воздействия в системе управления установкой выращивания монокристаллов.
1. Bagdasarov, X.S. Vysokotemperaturnaya kristallizaciya iz rasplava / X.S. Bagdasarov. - M.: Fizmatlit, 2004. - 160 s.
2. Sinel'nikov, B.M. Issledovanie kinetiki vyraschivaniya krupnyh monokristallov Al2O3 metodom Kiropulosa / B.M. Sinel'nikov, A.Yu. Ignatov, S.V. Moskalenko, O.T. Gudnin // Hi-miya tverdogo tela i sovremennye mikro- i nanotehnologii: VIII mezhdunar. konf. - Stavro-pol': SevKavGTU, 2008. - 458 s.
3. Es'kov, E.V. Opticheskiy sapfir, tehnicheskie trebovaniya i tehnologiya vyraschivaniya kri-stallov / E.V. Es'kov, A.Yu. Ignatov, V.S. Po-stolov, A.S. Filimonov // Himiya tverdogo tela i sovremennye mikro- i nanotehnologii: VIII mezhdunar. konf. - Stavropol': SevKavGTU, 2008. - 458 s.
4. Yudin, A.V. Osobennosti regulirovaniya moschnosti na vol'framovom vysokotemperaturnom nagrevatele / A.V. Yudin // Vestnik Rybinskoy gosudarstvennoy aviacionnoy tehnologiche-skoy akademii im. P.A. Solov'eva: sb. nauch. tr. - Rybinsk, 2011. - № 2 (21). - S. 108-112.
5. Lobacevich, K.L. Bespoiskovaya adaptivnaya sistema stabilizacii skorosti kristallizacii monokristallov / K.L. Lobacevich, A.V. Yudin // Avtomatizaciya i sovremennye tehnologii. - M.: Mashinostroenie, 2010. - № 11. - S. 23-26.
6. Gayduk, A.R. Algoritmicheskoe obespechenie adaptivnyh sistem upravleniya s identifikaciey / A.R. Gayduk, V.V. Shadrina // Vestnik Ivanovskogo gosudarstvennogo energeticheskogo universiteta. - 2018. - № 3. - S. 47-56.
7. Kuznecov, V.E. Adaptivnoe upravlenie s ekzomodel'yu tehnicheskim ob'ektom s ogranichennoy neopredelennost'yu / V.E. Kuznecov // Izvestiya SPbGETU «LETI». - 2016. - № 6. - S. 53-60.
8. Mel'nik, A.A. Postroenie modeley s primeneniem paketa System Identification Toolbox matrichnoy laboratorii MATLAB+SIMULINK / A.A. Mel'nik, E.A. Medunova // Sovremennye tehnologii. Sistemnyy analiz. Modelirovanie. - Irkutsk, 2015. - № 3 (47). - S. 93-97.
9. Fedosenkov, B.A. Teoriya avtomaticheskogo upravleniya: klassicheskie i sovremennye razdely: ucheb. posobie / B.A. Fedosenkov. - Kemerovo: KemGU, 2018. - 322 s.
10. Potemkin, V.G. Sistema MATLAB: sprav. posobie / V.G. Potemkin. - M.: Dialog - MIFI, 1997. - 350 s.
