Расчетно-экспериментальное исследование влияния эксплуатационного износа на динамические свойства кориолисового расходомера

Кориолисов расходомер (КР) предназначен для измерения массового расхода жидкостей и газов. Обзор конструкций КР приведен в [1–3], описание принципа их работы приведено в [4].

Представленный на рис. 1 КР состоит из основания 1, корпуса 2 и двух дельтаобразных трубок 3, по которым протекает рабочая среда. Под действием задающей катушки 4 трубки совершают изгибные колебания из плоскости ХY в противоположном друг к другу направлении на резонансной частоте, пропорциональной плотности рабочей среды. Поступательное движение рабочей среды во вращательном движении трубок приводит к возникновению кориолисовых сил, которые закручивают трубки вокруг оси Y и приводят к возникновению временной задержки Δ t между сигналами двух измери-

Рис. 1. Кориолисов расходомер: 1 – основание; 2 – корпус; 3 – трубки; 4 – задающая катушка;

5, 6 – измерительные катушки

тельных катушек 3 [4]. Эта задержка прямо пропорциональна массовому расходу.

Погрешность измерения массового расхода и плотности серийных КР достигает 0,1–0,05 %. Однако, при использовании КР для измерения расхода агрессивных или абразивных рабочих сред может возникать износ, коррозия и эрозия трубок, что приводит к увеличению погрешности измерения расхода и плотности [5, 6]. Кроме того, накопление на стенках трубок отложений, плотность которых отличается от плотности рабочей среды, также негативно влияет на точность измерения плотности.

Чтобы оценить техническое состояние КР, обычно его демонтируют и поверяют, проливая через него жидкость на специальном стенде [7]. Это дорогостоящая процедура, поэтому все чаще ₅₀                                                               Bulletin of the South Ural State University

Ser. Mathematics. Mechanics. Physics, 2026, vol. 18, no. 3, pp. 50–59

Яушев А.А.,                              Расчетно-экспериментальное исследование влияния

Логиновский В.А.                      эксплуатационного износа на динамические свойства … прибегают к методам поверки на месте эксплуатации, используя встроенные в электронику КР алгоритмы. Эти методы основаны на анализе изменения жесткости системы в процессе эксплуатации КР, что позволяет выявить возможный износ трубок.

Согласно [8], параметр жесткости k можно выразить через калибровочный коэффициент C F для определения массового расхода

„   rm    масса   ( сила/ускорение )       сила ₁

Cf =   =? = ----------7--- =------------= k •

A t ( время ) 2       ( время ) 2        перемещение

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

а = где mтрубки – масса трубки; Vтрубки – внутренний объем трубки, заполненный рабочей средой плот- ностью ρсреды

Р среды

k

т < mтрубки

^а 'трубки

трубки

В статьях [8, 9] и патенте [10] для оценки изменения жесткости предлагают: записать сигнал с измерительной катушки и сигнал возбуждения, и получить значения частотной передаточной функции FRF ( Frequency Response Function ) на резонансе и нескольких околорезонансных частотах; подобрать параметры модели с одной степенью свободы, чтобы она описывала основной режим работы колебательной системы КР; записать начальную жесткость модели. Для того чтобы удостовериться, что КР находится в пределах калибровки, полученную при каждой поверке жесткость сравнивают с начальной жесткостью.

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

Одним из существенных недостатков описанных методов поверки является сильная зависимость жесткости от температуры [13, 14], что повышает требования к точности калибровки температурных коэффициентов и может привести к ошибочным результатам поверки, если температура определена неправильно. Например, это может произойти при неравномерном распределении температуры вдоль измерительных трубок и в области установки датчика температуры.

В статье [15] и патенте [16] для оценки изменения жесткости предлагают измерять две собственные частоты КР и вычислять их отношение, которое представляет собой коэффициент, описывающий исходное состояние КР. Далее в процессе эксплуатации для оценки текущего состояния КР повторно измеряют две собственные частоты, сравнивают их отношение с исходным коэффициентом и делают вывод о состоянии системы. Основным преимуществом такого метода является существенно меньшее влияние температуры на результат поверки. На примере консольной балки показано, что отношение частот не зависит от модуля Юнга, следовательно, не зависит от его температурных изменений

2 ЕJ

P i = A²

P 2 A 2 ’

p i    ^ i      /4 , i ^1, 2 ^, 3...

ql ⁴

где λ – корни частотного уравнения; Е – модуль Юнга; J – момент инерции сечения; q –погонная масса балки; l – длина балки.

Для более сложной системы, представляющей собой консольную балку, разбитую на участки разного сечения (имитация локального износа) и с массой на конце (имитация навесных элементов на трубках), в патенте [16] показано, что отношение частот зависит от отношения массы груза к массе балки, отклонения момента инерции сечения и погонной массы каждого участка от номинала, но не зависит от изменения модуля Юнга, так как при его изменении частоты меняются одинаково. Полученные выводы проверены в статье [15] путем моделирования износов отдельных участков измерительных трубок при разных температурах с помощью конечноэлементной модели КР. Но в статье [15] рассмотрена только одна плотность среды и ограниченный набор износов отдельных участков трубки. Изменение плотности среды приводит к изменению отношения массы навесных элементов к массе трубки, что приводит к изменению отношения частот, которое необходимо оценить и учесть.

Данная статья дополняет рассмотренный в статье [15] и патенте [16] метод диагностики технического состояния КР. Выполнен расчетный анализ эффективности применения этого метода диагностики состояния КР при моделировании износов разных участков трубок для разных плотностей сред при разных температурах. Приведены результаты экспериментальной поверки.

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

Отношение частот зависит от плотности рабочей среды ρсреды, представим квадрат отношения рабочей и кориолисовой частот в виде уравнения f 2 = C1 раб/кор Рсреды + C2раб/кор, fкор где С1раб/кор и С2раб/кор – калибровочные коэффициенты, определяемые по рабочей и кориолисовой частотам при различных плотностях среды. Коэффициент С1 показывает, как сильно отношение частот зависит от плотности среды, геометрии трубок, расположения навесных элементов и их массы. Плотность рабочей среды ρсреды определяется по рабочей частоте КР, методика определения приведена, например, в патенте [16].

Для исследования влияния износа трубок на погрешность измерения массового расхода и собственные частоты КР используется конечноэлементная модель, приведенная в статье [15]. Она основана на предложенном Stack C.P. [17] методе моделирования с использованием балочных конечных элементов в формулировке Тимошенко, которые определяются через локальные матрицы масс, жесткости, демпфирования и гироскопических сил, учитывающие влияние текущей среды (массу среды, возникающие кориолиовые и центробежные силы) [18].

[M]{x(t)} + ([K] + [H(v)]){x(t)} + [C] {x(t)} = {P(t)}, где [M], [C], [K] и [Н(v)] – глобальные матрицы масс, жесткости, демпфирования и гироскопических сил формируются на основе локальных матриц конечных элементов модели, {x(t)} и{P(t)} – вектор-столбцы перемещений и внешних сил, v – скорость текущей среды.

Конечно-элементная модель колебательной системы КР (рис. 2) включает в себя основание, жестко закрепленное посередине, и две трубки с точечными массами, моделирующими катушки и перемычки между трубками. Текущая жидкость считается несжимаемой, однородной и невязкой, скорость течения жидкости по сечению и длине трубки постоянна.

Каждая трубка модели разбита на десять участков, расположенных симметрично слева и справа. Схема расположения участков на первой трубке приведена на рис. 3. Для каждого участ- ка имеется возможность независимо изменять диаметры, что позволяет проводить численные эксперименты и моделировать влияние основных факторов (износ трубок, в том числе локальный, отложения и др.), на массовый расход и собственные частоты КР. Влияние температуры в модели учтено следующим образом: в модель введена зависимость модуля Юнга от температуры; при изменении температуры учтено изменение линейных размеров поперечного сечения, обусловленного линейным расширением материала трубок, изменение внутреннего объема трубок и, соответственно, массы находящейся внутри трубок текущей среды.

Критерий k 1 диагностики состояния КР по отношению рабочей и кориолисовой частот выглядит следующим образом

^л 1    рраК/кор ^,

Л „ =ri „ n . +СР

^ раб/кор    C ¹ раб/кор ^p среды ' C ² раб/кор

22 раб кор ^.

Рис. 2. Схема расположения узлов колебательной Рис. 3. Схема разбивки первой трубки модели системы КР                                       на участки

Погрешность определения массового расхода

_{δ t =} Δ t ( T ) -Δ t ₀ ( T ) _{⋅ 100%,}

Δ t 0 ( T )

где Δ t – временная задержка при текущем состоянии КР, Δ t 0 – временная задержка при номинальных параметрах расходомера, T – температура.

Чтобы проверить эффективность предлагаемого метода диагностики состояния КР, рассмотрены следующие варианты износа: равномерный износ всех участков; износ отдельных участков трубки; симметричные износы участков, оказывающих максимальное и минимальное влияние на задержку; сочетания двух участков с максимальным и минимальным влиянием. Износы подобраны таким образом, чтобы они вызывали погрешность определения массового расхода примерно в 0,14 %. Погрешность определения массового расхода оценивается по изменению временной задержки относительно номинальных диаметров всех участков трубки. В табл. 1 приведены участки с износом и величина износа, выраженная как увеличение внутреннего диаметра трубки от номинала (8 мм). Для каждого варианта износа выполнены расчеты для трех сред при разной температуре (табл. 2).

Таблица 1

Положение и величина износа разных участков, выраженная через изменение внутреннего диаметра

Один участок		Симметричные участки		Два участка min & max		Все участки
T1_L1	8,04	T1_L3 & T1_R3	8,24 & 8,24	T1_L1 & T1_L3	8,022 & 8,182	All	8,004
T1_L2	8,22	T1_L3 & T2_L3	8,22 & 8,22	T1_L1 & T1_R3	8,022 & 8,182	–	–
T1_L3	8,33	T1_L1 & T1_R1	8,02 & 8,02	T1_L1 & T2_L3	8,022 & 8,182	–	–
T1_L4	8,144	T1_L1 & T2_L1	8,02 & 8,02	T1_L3 & T2_L4	8,04 & 8,14	–	–
T1_L5	8,058	–	–	–	–	–	–

Таблица 2

Плотности и температуры сред, используемые для расчетов

Среда	Воздух			Вода		Условное масло
Плотность, кг/м3	1,452	1,205	1,029	998,2	988,084	800	800	800
Температура, оС	–30	20	70	20	50	–20	20	90

Зависимость критерия k 1 диагностики технического состояния КР от температуры и плотности для различных вариантов износов трубки приведена на рис. 4. Для более детального анализа случаев износа, вызывающих небольшое изменение критерия, зависимость модуля критерия k 1 от температуры дополнительно приведена в логарифмическом масштабе на рис. 5.

Отклонение критерия k1 от нуля для номинального случая (отсутствия износов) в зависимости от температуры и плотности достигает порядка 10–6. Если в качестве референсного значения взять величину k1 = 10–5, то при использовании отношения двух частот (рабочей и кориолисовой) некоторые варианты износов окажутся недиагностируемыми. Поэтому дополнительно рассмотрена еще одна частота, соответствующая противофазной форме колебаний трубок (далее в тексте третья частота). Чувствительность каждого варианта отношения частот к износу отдельных участков трубки и равномерному износу двух трубок приведена на диаграмме рис. 6. Из нее видно, что отношение рабочей и кориолисовой частот для некоторых вариантов износа имеет избыточную чувствительность, повышающую вероятность ложного срабатывания сигнала о неисправности КР. Поэтому далее рассмотрены два критерия: k2, представляющий собой сумму модуля отклонения отношения рабочей и третей частоты и отношения кориолисовой и третьей ча- стоты и критерий k3, представляющий собой сумму модуля отклонения трех вариантов отношения частот. Результаты приведены на рис. 7–9.

& раб/III = C1 раб/III^рсреды + C ² раб/III - f раб /fill ;       & кор/III = C ¹ кор/III^рсреды + C ² кор/III - ^fKop / f ill •

7,0E-03 6,0E-03 5,OE-O3 ^ 4,0E-03 « & 3,0E-03 s. 2,0E-03 I.0E-03

k 2

| &раб/III

I + & кор/111 ;        k 3 =

& раб/кор | + | & раб/III | + < 1   па

кор/III

.

О ♦ — о ▲ •

>

<

О ♦ о A •

О ♦ о A •

о ♦ о ▲ •

•    1 1_LI                    /,VL-VJ •    T1L2              6,0Е-03 oTLL3              5,0Е-03 T1 L4 -               - 4,0Е-03 •T1 L5 S 3,0Е-03 -АН             о. OT1L3&T1R3 | 2,0Е-03 ♦ T1_L3 & T2_L3 И 1,0Е-03 ♦ T1 LI &Т1 R1      опг

> У

4 и

0,0E+00 -l,0E-03 -4

A—A

A

• 4

• _ф_

•

OT1L1&T2L1

8

o 5

>

О

О

о

AT1L1 &Т

L3        - 1 ,VC-UJ -R3       -2,0Е-03 _ТЗ L4

о -:

0

All LI & 1 20   40   60   80  100 iT1-L1&T. Температура, °C              T1 _L3 & T а)

200    400    600    800    1000 Плотность, кг/м3 б)

Рис. 4. Зависимость критерия k 1 диагностики технического состояния КР от температуры (а) и плотности (б) для различных вариантов износов трубок

Температура, °C

•Т1_Ы

• •    T1_L2 oTl_L3

T1_L4

• •    T1_L5

  -A11

OT1_L3&T1_R3 ♦ T1_L3 & T2_L3 ♦ T1_L1 &T1_RI OT1_L1 &T2_L1 ДТ1_Ы &T1.L3 aT1_L1 &T1_R3 ДТ1_Ы &T2_L3 a T1_L3 & T1_L4 ж Nominal

в браб/кор я ёраб/111 и § кор/Ш

Расположение участка с износом

Рис. 6. Влияние износов различных участков трубки на отношения частот

s

Рис. 5. Зависимость критерия k 1 диагностики технического состояния КР от температуры в логарифмическом масштабе для различных вариантов износов трубок

1,8E-O3

1,6E-O3 1,4E-03 L2E-03

8,0E-04

6,0E-04

4,0E-04

2,0E-04

0,0E+00

s

а)                                                      б)

Рис. 7. Зависимость критерия k 2 диагностики технического состояния КР от температуры (а) и плотности (б) для различных вариантов износов трубок

l,0E-02 1,0E-03 i l,0E-04 Qi О s a, *   1.0E-05 L0E-06 -4

*

ф

•    T1_L1                L0E-02 •    T1L2 OT1_L3 T1_L4                 l,0E-03 •    T1_L5 "Al1              ® OTIL3&T1R3   §. LOE-04 ♦    T1_L3&T2_L3    | ♦    T1_L1&T1_RI   §

* I * A

1

— 8 A

$ i •

8 *

11 • ♦

6 •

*

0 •

A ♦ <

►

A ♦

A ♦

▲ ♦

-

■

-

-

-

■ 0 -2

0

• * • 20   40   60   8 Температура, °C

о 1 1 LI Ж 14 LI         l.UL-UJ 4T1_L1 &T1_L3 aTILI &T1R3 *       Tl LI &T2 L3       L0E-06 0 100      “      “ a T1_L3 & T1_L4 • Nominal

• 0   -2

0

■ • • 20   40   60   8 Температура, °C

• 0   100

а)                                                       б)

Рис. 8. Зависимость критерия диагностики от температуры в логарифмическом масштабе для различных вариантов износов трубок: а – критерий k 2; б – критерий k 3

• Tl LI • T1_L2 о T1_L3

Плотность, кг/м3

TIL4 • T1_L5 -All

OTIL3&T1R3 ♦ T1_L3&T2_L3 ♦ TI_L1 &T1_R1 OTIL1 &T2L1 aTILI &T1L3 aTI LI &T1R3 AT1_L1 &T2_L3 ATI L3&T1 L4

а)                                                       б)

Рис. 9. Зависимость критерия k 3 диагностики технического состояния КР от температуры (а) и плотности (б) для различных вариантов износов трубок

На величину критериев k 1 , k 2 и k 3 оказывает влияние погрешность измерения плотности. Разность между эталонной плотностью и плотностью, определенной по рабочей частоте в зависимости от износа разных участков трубок, температуры и плотности среды, приведена на рис. 10.

• •    T1_L1

• •    Tl L2 oTl_L3

T1L4

• •    T1L5 -All

OTI L3 &T1_R3 ♦ T1L3&T2L3 ♦ T1_L1 &T1R1 OTI LI &T2JJ aTILI &T1L3 aT1_L1 &T1_R3

• •    Nominal

а)                                                       б)

Рис. 10. Разность между эталонной плотностью и плотностью, определенной по рабочей частоте в зависимости от износа разных участков трубок, температуры (а) и плотности (б) среды

Для проверки полученной зависимости погрешности измерения расхода и критериев диагностики технического состояния КР от локального изменения жесткости измерительной системы выполнены экспериментальные исследования. С помощью специального стенда (рис. 11) получен локальный электрохимический износ внутренней поверхности трубки (рис. 12). КР подсоединяли к оснастке, электрод помещали в измерительную систему к месту расположения перемычки, что примерно соответствует участку T1_L1. Включали насос, систему заполняли раствором NaCl 25 % и подавали ток с силой 4А в течение одной минуты. После этого КР относили на проливочную установку и определяли погрешность измерения расхода, которая составила примерно 0,03 %. Измеряли собственные частоты КР и определяли три критерия диагностики k 1 , k 2 и k 3 . Далее операцию электрохимического износа повторяли до достижения уровня погрешности примерно 0,07 %. Повторно измеряли собственные частоты КР и определяли три критерия.

Погрешность расхода, %

• к_1 Дк 2 хк 3

Рис. 13. Экспериментальная зависимость погрешности измерения расхода от трех критериев диагностики

Рис. 12. Фото внутренней поверхности разрезанной трубки после электрохимического износа

а)                                                      б)

Рис. 11. Специальный стенд для локального электрохимического износа внутренней поверхности трубок: а – фото стенда; б – схема стенда

На рис. 13 приведена экспериментальная зависимость погрешности измерения расхода и трех критериев диагностики для исходного состояния КР и после двух операций износа. Судя по величинам погрешности расхода и критериев, электрод (локализация места износа) находился между первым T1_L1 и вторым T1_L2 участками трубки (см. рис. 3).

В результате проведенного расчетно-экспериментального исследования применения отношения двух собственных частот КР для диагностики его технического состояния установлено:

• 1.    Использование отношения частот в качестве критерия диагностики позволяет обнаруживать недопустимое изменение погрешности измерения расхода, вызванное износом как локальных участков трубок, так и равномерным износом трубок.

• 2.    Износ различных участков трубки в разной степени влияет на погрешность измерения расхода и отношение частот. Достоверность диагностики можно повысить при использовании отношения трех частот (отношение различных пар, критерии k 1 , k 2 и k 3 ), но использование отношения рабочей и кориолисовой частоты более простое в реализации.

• 3.    Отношения частот слабо зависят от температуры, так как для каждого рассмотренного случая износа три критерия диагностики при изменении температуры находятся практически на

Яушев А.А., Расчетно-экспериментальное исследование влияния

Логиновский В.А. эксплуатационного износа на динамические свойства … одной прямой. В некоторых случаях при износе горизонтальных участков трубки наблюдается небольшой разброс, но он обусловлен больше влиянием плотности, а точнее, нарушением калибровки по плотности из-за износов.

• 1.    Wang, T. Coriolis Flowmeters: a Review of Developments Over the Past 20 Years, and an Assessment of the State of the Art and Likely Future Directions / T. Wang, R. Baker // Flow Measurement and Instrumentation. – 2014. – Т. 40. – С. 99–123.

• 2.    Даев, Ж.А. Эволюция кориолисовых расходомеров / Ж.А. Даев // Машиностроение: сетевой электронный научный журнал. – 2016. – Т. 4, № 3. – С. 33–39.

• 3. Anklin, M. Coriolis Mass Flowmeters: Overview of the Current State of the art and Latest Research / M. Anklin, W. Drahm, A. Rieder // Flow Measurement and Instrumentation. – 2006. – Vol. 17, Iss. 6. – С. 317–323.

• 4.    Raszillier, H. Coriolis-Effect in Mass Flow Metering / H. Raszillier, F. Durst // Archive of Applied Mechanics. – 1991. – Vol. 61, no. 3. – P. 192–214.

• 5.    Cunningham, T. Verification of Coriolis Flow Meter Calibration: Theory and Practice, Including Lab and Field Results / T. Cunningham // Can. Sch. Hydrocarb. Meas., 2015.

• 6.    Zhang, J. Early Detection of the Wear of Coriolis Flowmeters Through in situ Stiffness Diagnosis / J. Zhang, Y. Yan, J. Liu, E. Jukes // IEEE Transactions on Instrumentation and Measurement. – 2021. – Vol. 70. – P. 1–10. – Art no. 3520310.

• 7.    3124.0000.00-01 МП. Государственная система обеспечения единства измерений. Счетчики-расходомеры массовые. ЭЛМЕТРО-Фломак. Методика поверки. – 2021. – 36 с.

• 8.    Cunningham, T. An IN-SITU Verification Technology for Coriolis Flowmeters / T. Cunningham // Proceedings of the 7th International Symposium on Fluid Flow Measurement. – 2009. – С. 1–5.

• 9.    Rensing, M. Coriolis Flowmeter Verification via Embedded Modal Analysis / M. Rensing, T.J. Cunningham // Structural Dynamics, Vol. 3. Conference Proceedings of the Society for Experimental Mechanics Series. – New York: Springer New York, 2011. – С. 851–860.

• 10.    Патент RU 2628661. Вибрационный расходомер, а также способы и диагностика для поверки измерителя / Т.Ж. Каннингэм, Д.Д. Каполнек, М.Д. Ренсинг, К.Д. Ларсен; правообладатель МАЙКРО МОУШИН (США); заявл. 19.12.2014; зарег. 21.08.2017.

• 11.    Патент RU 2323419. Система и способ диагностики расходомера кориолиса / Ч.П. Стак; правообладатель МАЙКРО МОУШИН (США); заявл. 29.09.2003; зарег. 27.04.2008.

• 12.    Патент RU 2773633. Способ оценки состояния измерительной системы кориолисового расходомера / С.Е. Лобашев, М.Д. Шилин, Е.А. Горюнов, А.А. Яушев, В.А. Логиновский, А.В. Жестков; правообладатель Общество с ограниченной ответственностью «Элметро групп»; заявл. 07.12.2020; зарег. 06.06.2022.

• 13.    Costa, F.O. Modeling Temperature Effects on a Coriolis Mass Flowmeter / F.O. Costa, J.G. Pope, K.A. Gillis // Flow Measurement and Instrumentation. – 2020. – Vol. 76. – P. 101811.

• 14.    Патент RU 2803043. Способ оценки состояния кориолисова расходомера для его поверки и/или диагностики / А.А. Яушев, В.А. Логиновский, А.В. Жестков; правообладатель Общество с ограниченной ответственностью «ЭлМетро-Инжиниринг»; заявл. 11.05.2023; зарег. 05.09.2023.

• 15.    Яушев, А.А. Диагностика технического состояния кориолисова расходомера на основе оценки изменения отношения двух его собственных частот / А.А. Яушев, В.А. Логиновский // Законодательная и прикладная метрология. – 2025. – № 4(196). – С. 61–68.

• 16.    Патент US 8229695B2 Meter Electronics and Methods for Geometric Thermal Compensation in a Flow Meter / A.R. Pruysen, J.A.S. Maartensdijk, R.B. Garnett; правообладатель Micro Motion Inc. (США); заявл. 28.09.2006; зарег. 24.07.2012.

• 17.    Stack, C. A Finite Element for the Vibration Analysis of a Fluid-Conveying Timoshenko Beam / C. Stack, R. Garnett, G. Pawlas // 34th Structures, structural dynamics and materials conference. – 1993. – С. 1552.

• 18.    Тараненко, П.А. Влияние демпфирования на амплитуды и фазы колебаний системы с гироскопическими связями / П.А. Тараненко, Д.В. Телегин // Вестник Южно-Уральского государственного университета. Серия «Математическое моделирование и программирование». – 2025. – Т. 18, № 1. – С. 65–78.