Разработка компьютерной модели основного элемента агрегата дозирования топлива Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 21, № 2, 2016

Разработка компьютерной модели основного элемента агрегата дозирования топлива

И. Ш. НАСИБУЛЛАЕВ

Институт механики им. Р. Р. Мавлютова Уфимского научного центра РАН, Россия Контактный e-mail: sp.ishn@gmail.com

Исследуются режимы работы основного элемента топливной автоматики — сер-вопоршня дозирующей иглы. Построены математическая и компьютерная модели элемента. Методами полного факторного вычислительного эксперимента и ортогонального центрального композиционного планирования определены существенные факторы и найдены функции отклика. Получены зависимости параметров элемента от изменения температуры. Подготовлена модель сервопоршня в виде элемента вычислительного стенда.

Ключевые слова: агрегат дозирования топлива, сервопоршень, вычислительный эксперимент, полный факторный эксперимент, ортогональное центральное композиционное планирование, вычислительный стенд.

Введение

Проектирование динамических систем представляет собой сложный процесс, в результате которого должны быть выполнены конкретные требования по функционированию как системы в целом, так и отдельных ее элементов. Поскольку все технические системы состоят из связанного набора стандартных элементов (труб, гидродинамических сопротивлений, поршней, заслонок и др.), моделирование работы этих систем (течение жидкости, движение механических элементов) можно разбить на три этапа:

1-й этап — исследование отдельных элементов и построение целевых функций, определяющих интегральные параметры работы этих элементов (например, средний расход топлива) по управляющему сигналу (приложенному к элементу градиенту давления) и зависящим от температуры параметрам (материальным параметрам топлива, механическим элементам);

2-й этап — построение вычислительного стенда, объединяющего отдельные элементы;

3-й этап — построение полной модели системы для качественной и количественной оценки вычислительного стенда.

Отметим что, с одной стороны, модель, описывающая с достаточной степенью точности совместную работу хотя бы двух подвижных элементов, при использовании известных методов численного расчета получается достаточно громоздкой. С другой стороны,

© ИВТ СО РАН, 2016

разработчики всей системы автоматического управления (САУ) двигателем предъявляют следующее требование — модель исполнительной системы (агрегата топливопи-тания) должна функционировать в реальном масштабе времени. В связи с этим возникает вопрос о необходимости применения специальных методов для решения подобного класса задач.

При разработке вычислительного стенда техническое устройство (ТУ) представляется в виде набора связанных элементов. Каждый элемент имеет набор входных параметров (например, управляющий сигнал, температура, физические параметры элемента) и выходные данные (так называемые функции отклика), которые могут быть входными параметрами для следующего элемента. Сам элемент представляет собой правило перевода входных параметров в функции отклика в виде аппроксимаций. Свойства системы могут изменяться в ходе работы устройства (например, за счет нагрева элементов в результате трения), что может привести к смене режима работы устройства и аварийной ситуации. Для компенсации этих изменений необходимо менять параметры управляющего сигнала.

Существуют различные методы по определению чувствительности отклика системы к изменению параметров, основанные на аппроксимации отклика системы в параметрической области. Аппроксимация, основанная на факторном вычислительном эксперименте (ФВЭ) [1], позволяет построить приближенные модели с использованием небольшого количества вычислений, а результирующая полиномиальная функция — свести задачу полномасштабного моделирования математических моделей к набору простых зависимостей. При необходимости точность аппроксимации можно увеличить, повысив степень полинома или разбив интервал по одному или нескольким параметрам на подынтервалы. При разработке ТУ применение данного метода позволит провести поиск оптимальных параметров, что необходимо при исследовании систем большой размерности. Выделяя набор входных параметров методами ФВЭ, можно определить требуемые для разрабатываемого устройства выходные данные в виде функций отклика в форме полиномов от входных параметров системы. Работа с этими функциями позволяет:

• определить существенные факторы системы (воздействия на систему, меняющие динамические характеристики элементов);

• определить оптимальные значения факторов, обеспечивающие требуемые выходные режимы работы устройства путем поиска локальных экстремумов функции отклика;

• построить элементы исследовательского стенда для разработки ТУ, позволяющие изучать поведение системы в реальном времени (это не всегда возможно при моделировании полной математической задачи).

В данной работе в качестве первого элемента, для которого в рамках 1-го этапа будет проводиться исследование, выбран сервопоршень дозирующей иглы, осуществляющей подачу топлива в камеру сгорания газотурбинного двигателя. Этот элемент является основным для САУ частотой вращения газотурбинного двигателя, так как именно с его помощью производится стыковка электронной и гидромеханической частей системы управления. Принцип работы данного элемента подробно описан, например, в работе [2].

Целью настоящей работы является исследование режима работы сервопоршня. Для этого рассматривается чувствительность отклика элемента к изменению управляющего сигнала и параметров элемента методами ФВЭ, а также строится компьютерная модель сервопоршня для вычислительного стенда, который будет построен в рамках 2-го этапа.

1. Математическая постановка задачи

Рассматривается движение поршня длиной Ь, радиуса К и массой тр внутри трубы круглого сечения бесконечной длины, заполненной жидкостью (рис. 1). Система приводится в движение за счет перепада давления Ард(Ь), где Ар = (р\-р2) — амплитуда перепада давления; р\ и р2 — значения давлений слева и справа от поршня соответственно; д(Ь) — периодическая функция (-1 ^ д(Ь) ^ 1). В данной работе д(Ь) = со^(шЬ), ш = 2^/. Здесь £ — частота колебаний перепада давления. Скорость движения поршня имеет только одну компоненту ур вдоль оси Ох.

Уравнение движения поршня описывается вторым законом Ньютона:

где Я

V

нЫ;

V

тР~& = ^Ар - - Р(ьV), (1)

площадь поперечного сечения поршня; кр — коэффициент жесткости пружи-отклонение пружины от равновесного состояния. Сила сухого трения равна

Р(^р) = ^о sgn(wV), где ^о — модуль силы трения. Положение поршня хр описывается уравнением

(И

= Vр.

В начальный момент времени поршень покоится, т. е. 1>р(0) = 0, и находится в состоянии равновесия хр = 0. На покоящийся поршень действуют перепад давления (с начала процесса) и упругость пружины (со второй половины первого периода). Поршень начинает движение, когда сила, действующая на поршень, превысит силу трения покоя При движении на поршень действует постоянная по модулю сила трения скольжения /о. Сила трения скольжения Г0 всегда меньше силы трения покоя Для твердых тел отношение определяется коэффициентом kf = Р0/Р\, лежащим в диапазоне 0.6- 1.0 [3]. Отметим, что для параметров, когда поршень значительную часть периода покоится, коэффициент kf может оказывать влияние на динамику колебаний.

Для исследования системы запишем уравнения в безразмерном виде. В качестве характерного размера выберем радиус поршня К, а для характерного времени — обратную частоту управляющего сигнала ш-1. Размерные величины примут вид (безразмерные величины записаны с тильдами):

х

Кх,

'Р1

О X

Ь

5Р

Р2

Я \

V

к

р

Рис. 1. Геометрия задачи: поршень — сплошная линия; элемент трубы — штриховая линия

Опуская тильды, систему уравнений (1) и (2) запишем в безразмерном виде: ^ = Ьрд{г) - Ъкхр - Ь1 щп{ур), ^ = ур.

Здесь Ьр = Ар/(ЬррКш2) — амплитуда давления для поршня, где рр — плотность материала поршня; Ьк = и^/ш2, где = л/кр/тр — собственная частота колебаний пружины;

= Р0/(р0ЗрЬКш2) — сила сухого трения. Величины, характеризующие материальные параметры и геометрию системы, размерны.

2. Факторный эксперимент

Рассматриваемая задача зависит от нескольких факторов. Если провести расчет в определенных точках параметрического пространства, то можно построить аппроксимацию отклика системы в произвольной точке параметрического пространства.

Обозначим естественное значение г-го фактора как Xi, минимальное и максимальное значения — ж™1П и ж™ах соответственно, центральное значение — = (ж™ах + ж™1П)/2, интервал — xli = (ж™ах — ж™1П)/2. Кодированным значением i-го фактора Xi является линейное аффинное преобразование интервала [ж™1П, ж™ах] в интервал [—1, +1]:

/_

у _ Хi Хi

Лi = -1-•

Xi

Восстановление естественного значения фактора по кодированному: xi = xliXi + х°.

При выборе порядка аппроксимации необходимо учитывать, что более высокий порядок повышает точность аппроксимации, но требует большего количества расчетов (и вычислительных ресурсов). Для качественной оценки влияния факторов строится линейная аппроксимация. Если для отдельных факторов точность аппроксимации недостаточна, то порядок функции отклика увеличивается.

Рассмотрим линейную аппроксимацию отклика системы Y в виде полинома, зависящего от к факторов и их парных взаимодействий:

к к Y = Ьо + ^ biXi + ^ Ъ13ХгХ3 •

i=l i=l,i< j

Введем фиктивный фактор Х0 = 1 и обозначим парные взаимодействия XiXj ^ Х\,

к\

где I = к+1, • • • ,т; т = Ск = ^^-jy (количество парных взаимодействий при q = 2).

Функция отклика примет вид

к+т

Y = ^ biXi. =о

Коэффициенты bi в функции отклика Y определяются из условия минимума отклонений значений функции отклика от естественного значения у методом наименьших квадратов:

N

Ф(Ьо, •••, Ьк+т) = ^2(Уи — Yu) = min,

u=l

где N — количество измерений (моделирований), а функция Уи считается при тех же значениях факторов, что и естественное значение уи. Условие экстремума

дФ(Ъ0,...Л+ш1 = -2у- . . • ^ = 0, . = + т,

дЪ,

N / к+т \

: X] (Уи - X )

и=1 \ г=0 )

где Xги — значение г-го фактора (от 0 до к + га) в и-м опыте (от 1 до N).

Решая полученную систему линейных алгебраических уравнений, находим коэффициенты Ьг. Система имеет единственное решение, когда число опытов равно

N = 1 + к + га.

В общем случае, если учитываются взаимодействия более высоких порядков (до q2), то общее количество взаимодействий равно

92

га = £ С1.

1=2

Общее число опытов с максимальной кратностью взаимодействий = к равно

к к к N = 1 + к + ^ С«к = С°к + С1 + ^ С«к = ^ С«к = 2к,

д=2 д=2 д=0

что соответствует количеству опытов полного факторного эксперимента (ПФЭ) с двумя уровнями (функция отклика вычисляется для двух значений каждого фактора).

2.1. Полный факторный эксперимент 2к

Полный факторный эксперимент для к факторов на двух уровнях (ПФЭ 2к) позволяет определить существенные факторы модели и используется в качестве ядра для более эффективных (при меньшем количестве опытов) методов более высокого порядка. Найдем значения коэффициентов функции отклика методом наименьших квадратов. Систему N уравнений:

N / N-1 \

X ( Уи - X )

и=1 \ г=0 )

_ Уи - > ^ Ъгхги ) Х1и = 0, з = 0,...,Х - 1,

преобразуем к виду

^ уи хзи = ^ хгихзи, 3 = 0,...,Ы - 1.

N М-1 N

Ги^

и=1 г=0 и=1

Левая часть уравнений для ] = 0 (по определению Х0и = 1) равна

N N

у-^ги = XУи■

и=1 и=1

Вторая сумма правой части уравнений:

X,

N

и=1

N

Е 1

и=1 N

Е Х^и

и=1 N

Х-гиХ^и —

и=1

N

2 — N.

и

N

Е х*

и=1

0,

Ех

^ и=1

г — 3 — 0; 0, [г — 0,] — 0] или [г — 0,] г> 0, ]> 0, г — э; г — ]> 0.

0];

С учетом этих вычислений получим функцию отклика в виде

^ 1 N N

1

3 > 0.

у — ЬгХг, Хо — Хои = 1, Ьо — Ъ3 — Уихзи,

г=0 и=1 и=1

Количество расчетов, необходимых для построения функции отклика, равно Щ — 2 к.

Для построения полинома можно воспользоваться матрицей ПФЭ [1], которая содержит номер эксперимента (первый столбец); все возможные комбинации значений основных факторов Хг, включая фиктивный фактор Х0 — +1 (один столбец для фиктивного фактора и к столбцов для основных факторов); парные взаимодействия ХцХг2 (Ск столбцов); взаимодействия высших порядков (С% столбцов, q — кратность взаимодействий); столбцы со значениями отклика эксперимента у (количество столбцов соответствует числу исследуемых откликов). В качестве проверки в матрицу можно добавить следующие строки: ^ Хги — свойство симметрии, ^ Х2и — условие нормировки.

2.2. Ортогональное центральное композиционное планирование

Полный факторный эксперимент позволяет качественно определить значимость факторов и их взаимодействий и оценить характер зависимости функции отклика от факторов. Для проверки погрешности функции отклика можно провести ее сравнение с результатами расчета в наиболее удаленной точке параметрического пространства (в силу свойства ротатабельности в этой точке можно ожидать наибольшую невязку) — в центральной точке. Высокая погрешность означает, что линейной аппроксимации недостаточно и необходимо расширить полином квадратичными слагаемыми. Самый простой способ — ПФЭ вида 3 к — не является оптимальным, поскольку требует значительного количества опытов. Приведем более экономичный способ построения полинома второго порядка.

Рассмотрим полином ПФЭ 2к:

к к-1 к У — Ьо + ^ьгХг + ^^Ъгзхгхз + ...,

г=1 г=1 ]~>г

где Ьг — влияние факторов; Ъгз (г — ]) — парные взаимодействия; под многоточием подразумеваем взаимодействия более высоких порядков (часто они малы, и этими слагаемыми можно пренебречь), дополним его квадратичными слагаемыми:

к к-1 к к ¥ — Ьо + ^ ЬгХг + ^ ^ ЬгзХгХз + ■ ■ ■ + ^ ЬггX2,

г=1 г=1 ]~>г г=1

Ьц — квадратичные взаимодействия.

Для определения влияния квадратичных слагаемых необходимо добавить следующие опыты: все кодированные значения равны нулю, кроме одного, соответствующего квадратичному слагаемому и равному —а и +а. Всего таких опытов будет 2к. Добавляем один опыт в центральной точке (все кодированные значения равны нулю).

Полный факторный эксперимент составляет ядро из Щ = 2к опытов, и из него можно исключить малые взаимодействия (количество опытов это не уменьшит, но сократит количество вычисляемых коэффициентов полинома). Общее количество опытов ортогонального центрального композиционного планирования (ОЦКП) составляет N = Щ + 2к + 1.

Квадратичные слагаемые всегда неотрицательны, следовательно, нарушают симметрию относительно центральной точки плана. Среднее по опытам значение квадратичных факторов равно

!> = 1 УМ + 2а2

N ^ ги N

и=1

Величина плеча

а = — N0).

2

Условие симметрии можно получить центрированием квадратичных слагаемых X? — @:

к к-1 к к к г = Со + ^ ЪгХг + ^^ ЪгзХг Х3 + ••• + ^ Ьгг (X? — /3), Со = Ьо + /3 ^ Ьгг.

г=1 г=1 ]>г г=1 г=1

Величина плеча а определяется из условия ортогональности

N

N 1-

— /3)(Х2га — ¡3) = 0 ^ а =у1- — N0).

и=1 '

Найдем коэффициенты полинома. Составляя систему алгебраических уравнений методом наименьших квадратов аналогично ПФЭ, найдем

1 N к N / N

со = Уи, Ьо = со — Ьгг, Ь3 = ^ упХ3и ^ Х3,а.

и=1 г=1 и=1 / и=1

Здесь кратные и центрированные квадратичные факторы обозначены как независимые переменные.

Отличие в коэффициентах ОЦКП и ПФЭ связано с тем, что ОЦКП не является нормированным. Матрица планирования ОЦКП включает матрицу ПФЭ со существенными кратными взаимодействиями и столбцы с центрированными квадратичными слагаемыми. Также добавляются расчеты в звездных и центральной точках.

3. Результаты

Построим модель элемента, содержащего поршень и пружину. Фактором, приводящим поршень в движение, является перепад давления в виде Ард(Ь) = Ар ооз(ш1 + >^о), где Ар — амплитуда; ш — частота внешнего давления; (ро = -п/2 — начальный сдвиг фазы.

XI

Др(*

Х2

Х3, Х4, Х5

У1, Уз

У4

^5

X

Рис. 2. Схема элемента

Начальный сдвиг фазы в рассматриваемой задаче не влияет на динамику поршня, поскольку движение является симметричным и периодическим. Остаются два параметра: амплитуда х\ и частота х2 внешнего перепада давления. На схеме элемента (рис. 2) эти параметры представляют собой внешнее воздействие и не зависят от конструкции элемента. Фактором, влияющим на динамику элемента, является изменение температуры (во время работы поршень может нагреваться на 100 °С и более). Зависящие от температуры параметры — величина силы сухого трения х3 и упругость пружины х4. Величина силы сухого трения х3 пропорциональна силе реакции поршня и трубы, при нагреве поршень расширяется, что приводит к росту силы реакции и, следовательно, силы трения. С повышением температуры упругость пружины х4 уменьшается. В работе используется неявная зависимость физических параметров от температуры для большей универсальности (иначе потребуется явное указание температурной зависимости физических параметров в элементе). Последний из рассматриваемых параметров — отношение значения силы трения скольжения к силе трения покоя (kf) х5 [4] (х5 зависит от времени неподвижности поршня, приводящего к залипанию).

Диапазон значений параметров элемента определяется из физического смысла — величины не могут быть отрицательными или равными нулю. Максимальная величина силы сухого трения х3 определяется из условия bf < Ьр (внешнего перепада давления достаточно, чтобы вывести поршень из положения равновесия), т. е. < (Ар/Ь)Ур, где Ур — объем поршня. Отметим, что отклик системы не будет гармоническим, поскольку существует промежуток времени, в котором выполняется условие [2]

/о > ^ р сов(£) — к тах(|хр|), (3)

при котором поршень покоится при максимальном отклонении от состояния равновесия тах(|хр|). Чем больше промежуток времени, в котором выполняется это условие, тем выше будет отличие формы отклика от формы управляющего сигнала.

При выборе функций отклика необходимо учитывать, что функций должно быть достаточно для составления факторов, определяющих управляющий сигнал для последующего элемента. Дополнительные функции отклика определяют интегральные характеристики элемента за период (например, расход топлива), позволяют проверить гипотезы (например, неизменность периода колебаний) и провести дополнительное тестирование (например, по свойству симметрии задачи). Количество опытов определяется наиболее "сложной" функцией отклика (т. е. требующей большего количества коэффициентов полинома), а не количеством функций отклика (по одному набору опытов можно построить любое количество функций отклика).

В рассматриваемой задаче выбираем следующие функции отклика: путь поршня за период (характеризует расход топлива за период) у\; сдвиг фазы между колебаниями поршня и управляющего сигнала у2 (в единицах 2п); максимальное отклонение поршня (определяет амплитуду управляющего сигнала на следующем элементе) у3; период

колебаний поршня (определяет частоту управляющего сигнала на следующем элементе) у4; доля периода нахождения поршня в крайнем положении (например, в конструкциях с присоединенной к поршню запирающей иглой) у5.

Расчеты проводились с использованием схемы Рунге — Кутты четвертого порядка. Шаг по времени выбирался так, чтобы погрешность вычислений не превышала 0.01 % (2 • 105 подынтервалов на период). Количество периодов выбиралось по критерию установления периодического течения (относительное изменение функций отклика за два последних периода не более 10-5) и для различных параметров составляло 5-13 периодов.

Для определения существенных факторов проведен ПФЭ для выбранных к = 5 параметров в следующих интервалах: Ар = 2000± 200 Па, f = 40± 4 Гц, Ff = 7± 0.7 Н, кр = 2000 ± 200 Н/м, kf = 0.95 ± 0.05. Использовались следующие геометрические и физические параметры: диаметр трубы И = 2К = 10-2 м; длина поршня Ь = 10-2 м; плотность поршня (алюминий) рр = 2700 кг/м3 [5].

Функция отклика в общем виде (с учетом всех возможных взаимодействий) записывается как

У =Ьо + ЪгХг + № + Ьз Хз + ЪАХА + Ь5X5 + ЬеХ^ + Ь7ХгХз + Ъ8 ХЛ + №Х5+ + Ь\оХ2Хз + ЬцХ2Х4 + Ь\2Х2 Х5 + Ь\зХзХ4 + Ь\4 Х3Х5 + Ь15Х^Х5 + Ь\бХ\Х2Хзз + + ЬпХхХ2Х4 + Ь\8ХхХ2Х5 + Ьх^ХхХзХ^ + Ь20 ХхХзХ5 + Ь2\Х1 Х4Х5 + Ь22Х2ХзХ^ + + Ь2зХ2ХзХ5 + Ь24Х2Х4Х5 + Ь25ХзХ4Х5 + Ь26 Х\Х2ХзХ4 + Ь27 Х\Х2ХзХ5 + + Ь2нХ\Х2Х4Х5 + Ь29Х\ХзХ4Х5 + ЬзоХ2Хз Х4Х5 + Ьз\Х\Х2 ХзХ4Х5.

В табл. 1 представлены коэффициенты полинома Уг до взаимодействий второго порядка, поскольку влияние взаимодействий более высоких порядков незначительно. Функция отклика изменяется в диапазоне [1.54, 3.35], т.е. на 37% от среднего значения 2.45. С учетом взаимодействий второго порядка максимальная погрешность полинома, рассчитанная на двух уровнях (т. е. в точках, по которым строилась аппроксимация -1,1), составляет 1.96%, на пяти уровнях (включая промежуточные точки -1,-0.5, 0, 0.5,1) — 3.69 %. Существенное воздействие на функции отклика у\ оказывают упругость пружины (положительная связь, т. е. увеличение параметра повышает значение функции отклика), амплитуда перепада давления (положительная связь) и частота управляющего сигнала (отрицательная связь). Влияние силы трения и отношения величин силы трения покоя и трения скольжения незначительны. Наиболее существенные парные взаимодействия — ХгХ4 и Х2Х4.

На сдвиг фазы У2 влияют амплитуда перепада давления, сила трения и отношение сил трения. Наибольший вклад вносят парные взаимодействия ХгХз, ХгХ5, ХзХ5. В целом изменение У2 не более 1.5%.

Наибольшее влияние на долю периода в состоянии покоя У5 оказывают амплитуда перепада давления, частота управляющего сигнала и упругость пружины вследствие

Таблица 1. Коэффициенты полиномов Уг, У2 и У5 ПФЭ

Полином Ьо Ъг Ь2 Ьз Ь4 Ь5 Ьб Ъ7 Ьн Ь9 Ъго Ьп Ь\2 Ь\з Ь\4 Ь\5

Уг • 102 235.7 27.2 22.7 -3.6 -34.8 -2.3 2.5 0.6 -3.9 0.4 -0.2 -3.2 -0.1 0.4 -0.6 0.3

12 • 104 5745.3 -79.0 7.3 79.0 -3.6 42.2 -0.5 -9.7 0.3 -5.3 0.5 -0.7 0.3 -0.3 5.3 -0.2

П • 102 288.0 -3.1 -21.4 3.1 10.8 -2.5 -0.5 -0.2 0.3 0.2 0.5 1.1 0.1 -0.3 -0.2 -0.1

соотношения (3). Наиболее существенные парные взаимодействия — Ь8XiХ4 и biiХ2Х4. Суммарный вклад — порядка 15%.

В рассматриваемой параметрической области колебания поршня являются симметричными относительно состояния равновесия, следовательно, Y3 = Yi/4 — это условие выполняется с точностью ошибок округления. Форма сигнала существенно зависит от величины силы сухого трения — с ее ростом форма сигнала меняется от гармонического до ступеньки (рис. 3). Если форма сигнала для связанного элемента имеет значение, то необходимо учитывать функцию отклика Y5. Период колебаний поршня Y4 определяется (в пределах вычислительной погрешности 0.01 %) периодом управляющего сигнала. Влияние собственных колебаний поршня при наличии сухого трения отсутствует [2]. Время установления периодического движения и амплитуда отклика существенно зависят от параметров.

Параметр Х5 имеет наименьшее влияние на функции отклика и может быть исключен из дальнейшего рассмотрения. Отметим, что в техническом устройстве можно полностью исключить параметр Х5 за счет приложения высокочастотного малоамплитудного сигнала. Расчеты показывают, что такой сигнал практически не влияет на динамику поршня, но при этом предотвращает залипание поршня к поверхности канала. Взаимодействиями параметров с кратностью выше двух также можно пренебречь.

Рассмотрим характер нелинейности искомых функций. Для этого построим матрицу ОЦКП, включающую матрицу планирования ПФЭ до взаимодействий второго порядка по четырем параметрам (к = 4, параметр Х5 не рассматривается) и центрированные квадратичные слагаемые Хц — ß. Тогда N0 = 24 = 16, N = N0 + 2к + 1 = 25, плечо a = 21/2, ß = 0.8 и полином примет вид

Y = со + bi Xi + b2 Х2 + b3 X3 + 64 X4 + b5 Xi X2 + b6 Xi X3 + b7 Xi X4 + b8 X2 X3+ + ь9 X2 X4 + bio x3 X4 + bii (X2 — ß) + bi2 (XI — ß) + bi3 (X32 — ß) + bi4 (X2 — ß).

(4)

Рис. 3. Зависимость координаты поршня хр от безразмерного времени tf (Х5 = 1): сплошная линия — Ар = 1800 Па, / = 36 Гц, = 6.3 Н, кр = 1800 Н/м; штриховая линия — Ар = 2200 Па, / = 36 Гц, ^о = 7.7 Н, кр = 1800 Н/м; пунктирная линия — управляющий сигнал

Таблица 2. Коэффициенты полиномов У г, У2 и У5 ОЦКП для больших значений ^о

Полином со Ьг Ь2 Ьз Ь4 Ь5 Ьб Ьг Ь8 Ь9 &10 Ьп &12 Ьгз Ьг4

п • 102 2327.5 274.8 225.5 -41.3 -346.3 24.8 7.1 -39.8 -2.1 -32.1 5.2 -3.8 -2.0 -3.7 42.9

У2 • 104 5785.2 -84.4 7.5 83.9 -3.8 -0.6 -10.5 0.3 0.6 -0.7 -0.3 9.6 0.5 1.0 0.4

П • 102 285.7 -2.9 -21.3 2.9 10.6 -0.5 -0.2 0.2 0.5 1.1 -0.2 0.2 0.0 0.0 -0.8

В табл. 2 представлены коэффициенты полиномов с учетом парных и квадратичных взаимодействий. Функция уг является нелинейной с наибольшим квадратичным вкладом от X\.

Максимальная погрешность полинома, рассчитанная на пяти уровнях (включая промежуточные точки -1, -0.5, 0, 0.5, 1), составляет 0.46%, что значительно ниже погрешности полинома ПФЭ. Для функции Уг существенным нелинейным параметром является упругость пружины. Влияние остальных нелинейных параметров на порядок ниже. Нелинейность функции У2 определяется амплитудой перепада давления, функция У5 преимущественно линейная.

Отметим, что при малых значениях силы трения (^0 < 1 Н) отклик системы остается гармоническим (на рис. 4 перемещение поршня при малой силе трения показано сплошной линией, при высоком значении силы трения — штриховой), т.е. форма сигнала похожа на форму управляющего сигнала. В табл. 3 представлены коэффициенты полинома в том же параметрическом пространстве, но для малой величины силы трения (^о = 0.5 -0.75). Видно, что форма отклика слабо отличается от управляющего сигнала (функция У5 мала). Сдвиг фазы практически не зависит от параметров и составляет половину периода (2пУ2 ~ ^). Наибольшее влияние на величину пути поршня за период Уг оказывают амплитуда перепада давления и жесткость пружины. Также стоит отметить,

| ц / ...»—|..г.. ;! I: V.......... 1 ■ й ; "ч *• А : : 1/ || : Г.; ! 1 \ \\! ! 1 I1, \ 7 || и V 1 \ Л ! 1 1 \ ■ ■ 7 Р • У п : % / 1 ! \ А : г 1 ; 11 I'«

и . ........;.л ■ 1 ........... :| •И \ ' и II 1.111.....ц "Лт ■ : V : \ 1 > \ > 1 1 ____.......Л"....... 1 1 |/ ', 1 ||1 /| 1 1 —г ; V 1 1 \ У \ • 1 ........V....... || 11 //;Ц Г ' '1 Ч • \ У \ ! .......л/......

0 1 2 3 4 5 6

Рис. 4. Зависимость координаты поршня хр от безразмерного времени £ • / (Х5 = 1) для Ар = 2000 Па, f = 40 Гц, кр = 2000 Н/м: оплошная линия — = 0.5 Н; штриховая линия — = 3 Н; пунктирная линия — управляющий сигнал д (£)

Таблица 3. Коэффициенты полиномов Y\, и ОЦКП для малых значений Fq

Полином Со bi Ь2 Ьз Ь4 b5 b6 Ьт b8 bg bio bii bi2 bi3 bi4

Yi ■102 3374.0 342.0 37.0 —9.6 —356.4 7.7 0.2 —38.8 —7.4 19.0 5.8 — 1.2 —23.0 —1.3 33.7

12■104 4988.3 1.2 0.4 1.3 8.8 —3.2 —2.8 —0.3 1.6 —0.9 — 1.3 —3.7 —9.8 —3.3 —5.3

Y5 ■ 102 7.0 — 1.2 —8.2 0.2 —3.7 1.7 0.7 0.1 —0.4 5.1 0.3 3.2 4.6 4.2 3.2

что периодический режим движения поршня устанавливается за большее количество периодов, чем при движении с большей величиной силы трения. Максимальная на пяти уровнях погрешность полинома составляет 0.75%.

Погрешность аппроксимации сильно зависит от диапазона варьирования параметров. Например, увеличение диапазона варьирования с ±10 до ±15 % в рассматриваемой задаче приводит к росту погрешности с 0.5 до 15 %. Точность полинома можно повысить несколькими способами. Во-первых, увеличить количество уровней варьирования за счет нелинейного увеличения числа опытов. Во-вторых, можно повысить порядок полинома, если источником погрешности является нелинейное поведение функции отклика. В-третьих, можно разбить интервал варьирования на несколько подынтервалов и построить функции отклика для каждого. Количество опытов в этом случае увеличится пропорционально числу интервалов. Последний вариант удобно использовать при расширении интервала по одному параметру. При этом получаются две функции отклика для каждого интервала. При аппроксимации значение функции в точке не совпадает со значением полинома в этой точке и функция отклика Y, полученная для двух интервалов, может оказаться разрывной на границе интервала. Этого можно избежать, если взять частично пересекающиеся интервалы и на их пересечении (Ха, Хь) построить линейную интерполяционную сшивающую функцию

X — X

Y(X) = ВД) + [Yl(Xa] - ВД)].

Если текущий (базовый) элемент связан со следующим элементом, то необходимо построить соответствие между функциями отклика и входными параметрами связанного элемента. Частота управляющего сигнала не меняется, следовательно, х2 базового элемента соответствует х2 следующего элемента. Амплитуда колебаний, форма сигнала и сдвиг фазы меняются. В базовом элементе управляющий сигнал имеет вид

Ap(t) = Yi cos [2n(x2t + ро)].

В связанном элементе для малой величины силы трения (до F0 = 1 Н)

Ap(t) = хз cos [2n(x2t + ро + ^2)] •

Для больших значений сил трения (выше F0 = 1 Н) необходимо построить симметричную кусочно-линейную ступенчатую функцию:

• Ap(t) = — Хз на интервале 2п[у2 — ^о; У2 — фо + Уь/2];

• Ap(t) = +х3 на интервале 2ж[у2 — уо + 1/2; У2 — + 1/2 + Уъ/2];

• в остальной области интервалы соединяются линейной функцией.

Проведем оценку диапазона изменений параметров элемента при изменении температуры. При движении поршня в результате действия силы трения F0 происходит

нагрев контактирующих поверхностей. За один период колебаний (г) поршень проходит путь, равный ух, при этом совершается работа силы трения А = Ух Я0, переходящая в тепло Я = А/2 (половина тепла переходит в поршень, половина — в стенку трубы). Среднее изменение температуры поршня равно

8Т Я ВД

трс

где тр = 2.1 г — масса поршня; о€ = 0.92 кДж/(кг-К) — удельная теплоемкость. Для среднего значения ух = 0.024 м и силы трения Я0 = 1 Н получим изменение температуры за один период и 1 с:

5Т(т) = 0.00615 К, 5Т(1) = 5Т(т)/ = 0.25 К.

При нагреве поршня происходит относительное удлинение радиуса [6]:

^ = А«Г, (5)

где Л — коэффициент линейного теплового расширения. Стенка трубы компенсирует это расширение за счет напряжений. На поверхности возникает сила реакции, равная

= —ЕаЗо, К

где Еа = 7 • 1010 Па — модуль Юнга для алюминия; — эффективная площадь контакта, равная средней за период площади контакта поверхностей с разными температурами. В начальный момент времени эта площадь равна нулю (поршень и труба имеют одинаковую температуру). При перемещении поршня на ¿х эффективная площадь увеличивается на 2пКс1х. За период колебаний в0 ~ 2пЯЬУ1/(У1 + V). С учетом (5) и закона Амантона — Кулона (сила сухого трения пропорциональна силе реакции опоры Дх3 = ^^р, где ^ = 0.4 — коэффициент трения) получим поправку к силе трения

Дхз = /лХЕавобТ(г) « 0.9 Н. (6)

Эта оценка является оценкой сверху, поскольку не учитывает явление теплопереноса. Величина силы сухого трения с учетом изменения температуры равна х3 + Дх3.

Коэффициент жесткости пружины кр = Е3в/Ь3, где Е3 — модуль Юнга для стали; в и Ь3 — площадь сечения и длина пружины. Геометрия пружины зависит от температуры незначительно (изменение в/Ь3 пропорционально коэффициенту линейного теплового расширения и имеет порядок 10-5). В температурном интервале 20-160 °С модуль Юнга линейно уменьшается на 7.5% [7], следовательно, зависимость от температуры параметра х4 можно записать в виде

х4(Т) = х4(Т = 20 °С)0.075[1 + (Т - 20 °С)/140 °С]. (7)

Таким образом, выбранные интервалы параметров являются адекватными.

Рассмотрим пример использования компьютерной модели. Сначала определяются коэффициенты функции отклика (4) и существенные параметры и взаимодействия. Рассчитывается зависимость температуры от времени, и обновляются значения параметров (6) и (7). Для обновленных параметров с использованием полинома вычисляются функции отклика (4). По функциям отклика строится управляющий сигнал для связанного элемента.

Заключение

Целью данной работы является построение элемента вычислительного стенда для основного элемента топливной автоматики — сервопоршня дозирующей иглы, управляемого периодическим перепадом давления, с учетом упругости пружины и силы трения.

Решение задачи механики осуществляется методами численного моделирования неоднородных обыкновенных дифференциальных уравнений (схема Рунге — Кутты 4-го порядка). Анализ результатов компьютерного моделирования проведен методами факторного вычислительного эксперимента (полного факторного эксперимента и ортогонального центрального композиционного плана).

Построены математическая и компьютерная модели сервопоршня. Выбраны факторы вычислительного эксперимента и функции отклика. Предложена аппроксимация функций отклика в виде полинома с линейными, парными и квадратичными вкладами факторов. Определены существенные факторы динамики поршня и их взаимодействия и погрешность аппроксимации. Установлены температурные зависимости параметров, что позволяет провести адекватное моделирование работы элемента с учетом динамических изменений в этом элементе в процессе работы. Показан способ объединения отдельных элементов вычислительного стенда в систему.

Показано, что за конечное количество периодов перепада давления (от 5 до 13) устанавливается периодический и симметричный отклик системы (перемещение поршня). Форма отклика в области больших давлений имеет вид ступеньки, с уменьшением величины силы трения скольжения форма становится более гармонической (близкой к форме перепада давления). На динамику движения поршня наибольшее влияние оказывают упругость пружины, амплитуда перепада давления с положительной связью и частота перепада давления с отрицательной связью. На сдвиг фазы перемещения поршня относительно перепада давления влияют амплитуда перепада давления, сила трения и отношение сил трения.

Представлены: способ построения управляющего сигнала для связанного элемента стенда по функциям отклика; схема расчета скорости нагрева из-за работы силы трения (порядка 0.25 град. в секунду); дополнительная величина силы трения (порядка 1 Н), возникающая за счет теплового расширения; влияние нагрева на изменение жесткости пружины.

В работе предложена схема элемента вычислительного стенда сервопоршня, связывающая управляющий сигнал (форма внешнего перепада давления) с функциями отклика (динамика перемещения поршня) и учитывающая изменение параметров под действием температуры (изменение силы трения поршня и упругости пружины в результате нагрева) в ходе эксплуатации элемента. Главным преимуществом данного подхода является замена полного моделирования (требующего больших вычислительных ресурсов) аппроксимацией, позволяющей проводить расчеты функций отклика в реальном времени.

Предложенная схема построения компьютерной модели элемента универсальна и может быть использована для разработки моделей других элементов, а также для построения вычислительного стенда системы элементов.

Благодарности. Работа выполнена при финансовой поддержке РФФИ (гранты № 14-01-97019-р_поволжье_а, № 14-08-97027-р_поволжье_а).

Список литературы / References

[1] Адлер Ю.П., Маркова Е.В., Грановский Ю.В. Планирование эксперимента при выборе оптимальных условий. М.: Наука, 1976. 280 с.

Adler, U.P., Markova, E.V., Granovskiy, U.V. Experiment plans with selecting the optimal conditions. Moscow: Nauka, 1976. 280 p. (In Russ.)

[2] Денисова Е.В., Насибуллаева Э.ШШ., Насибуллаев И.ШШ. Исследование динамических процессов в элементах топливной автоматики // Мехатроника, автоматизация, управление. 2014. № 5. С. 31-36.

Denisova, E.V., Nasibullaeva, E.Sh., Nasibullayev, I.Sh. Investigation method of dynamics processes in the fuel elements of automation // Mekhatronika, Avtomatizatsiya, Upravlenie. 2014. No. 5. P. 31-36. (In Russ.)

[3] Современная трибология: итоги и перспективы / Отв. ред. К.В. Фролов. М.: ЛКИ, 2008. 480 с.

Modern tribology: results and prospects / Otv. red. K.V. Frolov. Moscow: LKI, 2008. 480 p. (In Russ.)

[4] Насибуллаев И.ШШ., Насибуллаева Э.ШШ., Денисова Е.В. Влияние различных видов силы трения в системе двух коаксиальных цилиндров // Мехатроника, автоматизация, управление. 2014. № 10. С. 54-59.

Nasibullayev, I.Sh., Nasibullaeva, E.Sh., Denisova, E.V. Influence of different types of friction forces in the two coaxial cylinders // Mekhatronika, Avtomatizatsiya, Upravlenie. 2014. No. 10. P. 54-59. (In Russ.)

[5] Алюминий: свойства и физическое металловедение: Справ. изд.; Пер. с англ. / Под ред. Дж.Е. Хэтча. М.: Металлургия, 1989. 422 с.

Aluminum: properties and physicsl metallurgy: American society for metals / Ed. by J.E. Hatch. Ohlo: Metals Park, 1984.

[6] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. VII: Теория упругости. М.: Наука, 1987. 248 с.

Landau, L.D., Lifshitz, E.M. Theory of elasticity. Vol. 7 (3rd ed.). Butterworth-Heinemann, 1986.

[7] ГОСТ 1050-88. Прокат сортовой, калиброванный, со специальной отделкой поверхности из углеродистой качественной конструкционной стали. Общие технические условия. Введ. 1991-01-01. М.: Стандартинформ, 2010. 19 с.

State Standard 1050-88. Rolled metal, calibrated, with special surface finish made from quality carbon structural steel. General specifications. Moscow: Standartinform Publ., 2010. 19 p. (In Russ.)

Поступила в 'редакцию 26 октября 2015 г., с доработки —15 декабря 2015 г.

The development of a computer model for the main element of the fuel metering unit

Nasibullayev, Ildar Sh.

Mavlutov Institute of Mechanics RAS, Ufa, 450054, Russia Corresponding author: Nasibullayev, Ildar Sh., e-mail: sp.ishn@gmail.com

The purpose of this paper is to construct an element of the computational stand for a servo-piston induced by periodic pressure gradient. The research takes into account the elasticity of the spring and friction forces.

© ICT SB RAS, 2016

Solution of the problem is done by numerical modelling of inhomogeneous ODE (Runge-Kutta method). Analysis of the results of computer simulation carried out by methods of computational experiment. Response functions and leading factors were found.

It was shown that after finite number of periods periodic and symmetrical response of the system is established. Leading factors of piston motion are elasticity of the spring, the amplitude of the pressure gradient and the frequency of the pressure gradient.

In the framework of computational stand we present a method of constructing a control signal for the associated element of the response function, a scheme for calculating the rate of heating due to friction, the correction of the frictional force due to thermal expansion as well as the change of the spring elasticity due to heating.

In this paper we introduced an element of a computer stand servo-piston connecting external factors (form of external pressure gradient) with the response function (motion of the piston) and changes in internal factors during the process. The main advantage of this approach is the replacement of the complete simulation by approximation, which allows calculations in real time.

The proposed scheme for constructing a computer model of the element is universal and can be used to develop models of other elements, and also to build the computational stand for the system of elements.

Keywords: fuel metering unit, servo-piston, computational experiment, a full factorial experiment, the orthogonal central compositional planning, computer stand.

Acknowledgements. The reported study was funded by RFBR (projects № 14-0197019 and № 14-08-97027).

Received 26 October 2015 Received in revised form 15 December 2015