OTIP

~(k )

pt

или

~

~

(H

(k )

)

3

∫∫

∂ϖm ∂W

sin θ(H

(k )

)

3

+

∂ϖm ∂W

(k )

=

∂θ ∂θ

∂ϕ ∂ϕ

sin θ

dϕdθ P

Ω( k )

	~
	∂ϖ		Λsin2
					( m =1, N × M ).			(10.23)
= − ∫∫		m		θ(ε1(k ) cosϕ + ε(2k ) sin ϕ) dϕdθ
Ω( k )	∂ϕ

На границе области k			функция давления задана краевыми условиями (10.13).
Неизвестными являются			значения функции в узлах, расположенных внутри
области Ω(k). Следовательно, искомые узловые значения функции давления в
области k содержатся в уравнениях (10.23) с номерами m =
				M(d −1) +1, Md −1;
d =		. Таким образом, получена система (N – 2)(M – 2) линейных алгебраи-
	2, N −1

ческих уравнений относительно	~(k )		, которая в матричной форме имеет вид
	pt
A(k )		P(k ) = C(k ) .		(10.24)

Компоненты квадратной матрицы A(k ) и матрицы-столбца С(k ) определяются по формулам:

Am,(k )t

Cm(k )

			~	~	~	~	(H	(k )	)	3	dϕdθ,
=	∫∫		∂ϖm	∂ϖt (H	(k ) )3 sin θ + ∂ϖm ∂ϖt
			∂θ	∂θ	∂ϕ	∂ϕ	sin θ
	Ω( k )
			~
= − ∫∫			∂ϖ	m Λsin2
					θ(ε1(k ) cosϕ + ε2(k ) sin ϕ) dϕdθ.
	Ω( k ) ∂ϕ

(10.25)

(10.26)

Так как весовые и пробные функции задаются на каждом конечном элементе формулами (10.16), то при нахождении двумерных интегралов (10.25), (10.26)

коэффициентов Cm(k ) , Am,(k )t внутри вычислительной области Ω(k) ненулевые вклады

вносят только элементы A, B, C, D, примыкающие к узлу m, совпадающему с узлом t или (i j) (см. рис. 10.7).

В результате дискретизации граничных условий (10.7) получаем

~(k )	~(k )
p j	=1, p(N −1)M + j =1 ( j =1, M );
~I	~II	~II			~I
p(i−1)M +1 = p(i−1)M +M , p(i−1)M +1				= p(i−1)M +M

(10.27)

(i = 2, N −1).

Две системы уравнений (10.24) при k = I и k = II и равенства (10.27) совместно представляют линейную алгебраическую систему уравнений для вычисления искомых сеточных значений функции давления в зазоре гидродинамического подвеса. В матричной форме полученную систему уравнений представим в виде

A P = C,

(10.28)

где P – столбец, составленный из искомых узловых значений функции давления = pi(kj ) ; C – столбец, состоящий из коэффициентов, обусловленных

дискретизацией краевой задачи; A – квадратная разреженная матрица почти девятидиагональной структуры с локальными отклонениями, элементы которой

171

обусловлены дискретизацией. Ненулевые элементы as q матрицы А расположены по указанным на рис. 10.8. диагоналям и в отмеченных (×) узлах.

	Решив		систему			линейных
	алгебраических			уравнений			(10.28),
	найдем	распределение				давления в
	зазоре	гидродинамического					подвеса
	шарового гироскопа.
	Для	решения			систем линейных
	алгебраических			уравнений			высокого
	порядка		с ленточной			структурой
	матрицы		А эффективным				является
	построенный			на	основе		метода
	Гаусса		обобщенный				алгоритм
	Томаса. Основная идея состоит в
	многократном повторении прогонки
Рис. 10.8. Структура матрицы А	вперед с		целью		последовательного
	исключения диагональных элементов

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

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

Результирующая гидродинамическая реакция подвеса F , как и в схеме А, равна

F = F I + F II .

Проекции Fxi (i =1, 3) результирующей гидродинамической силы на оси СКxi статора представляем в виде

Fxi = KF F xi ,

где KF = πpo R22 – масштабный коэффициент для сил; F xi безразмерные проекции реакции жидкости:

F x1 = − 1 π∫[π∫ pI cosϕ

π 0 0

F x2 = − 1 π∫[π∫ pI sin ϕ

π 0 0

F x3 = − 1 π∫[π∫ pI dϕ +

π 0 0

dϕ +

dϕ +

2π II

∫ p

π

	2π				II	cosϕ dϕ]sin2
			p
	∫						θ dθ,
	π
	2π				II	sin ϕ dϕ]sin2
	∫ p						θ dθ,
	π
	dϕ]cosθsin θ dθ.

= F Ixi + F IIxi –

(10.29)

172

Главный момент реакций жидкости M O , приложенный к ротору, равен сумме моментов M OI и M OII . Введя масштабный коэффициент для моментов

KМ = πμω R24 δ, записываем

MOxi = KM M Oxi = KM (M OxiI + M OxiII ) (i =1, 3).

Безразмерные величины возмущающих моментов для гироскопа MOx1 , и момента сопротивления его быстрому собственному вращению находятся в виде интегралов:

M Ox2 ,

MOx3

M Ox1

M Ox2

M Ox3

3

π π

∂

I

∂

I

H I

Λ sin θ

p

H I sin

p

=

{∫∫

ϕ+(

+

H I

)cosθcosϕ sin θ dθ dϕ+

Λπ

∂θ

∂ϕ

sin θ

3

0 0

II

II H II

2π π

∂

∂

Λ sin

θ

p

H II sin

p

+

∫ ∫

ϕ+(

+

H II

)cosθcosϕ sin θ dθ dϕ},

∂θ

∂ϕ

sin θ

3

π 0

3

ππ

∂

I

∂

I

H I

Λsin θ

p

H I

p

=

{∫∫

−

cosϕ+(

+

3H I

)cosθsin ϕ sin θ dθ dϕ+

Λπ

∂θ

∂ϕ sin θ

0 0

2π π

∂

II

∂

II

H II

Λ sin θ

p

H II

p

+

∫

∫ −

cosϕ+(

+

3 H II

)cosθsin ϕ sin θ dθ dϕ}, (10.30)

∂θ

∂ϕ

sin θ

π 0

3

ππ

∂

I

H I

Λ sin θ

p

=

{∫∫

−

−

H I

sin2 θ dθ dϕ+

Λπ

∂ϕ

sin θ

3

0 0

II

2π π

∂

H II

Λ sin θ

p

+

∫ ∫

−

−

H II

sin2 θ dθ dϕ}.

∂ϕ

sin θ

3

π 0

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

10.4. Алгоритмы определения характеристик гидроподвесов, выполненных по схемам А, В

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

173

1. Ввод параметров подвеса: μ – динамического коэффициента вязкости, pо – давления в камере, ω – угловой скорости вращения ротора; геометрических характеристик: R2 – радиуса ротора, δ – разности радиусов полости статора и ротора, величин относительных смещений центра ротора εi (i =1, 3 ) вдоль координатных осей СКxi; величины χ сдвига центров сегментов I и II; угла ориентации вектора сдвига φ; усечений полусфер I и II, составляющих статор, χiI , χiII (i = 3 – для схемы А , i = 2 – для схемы В); количества узлов сетки в областях I, II по переменным θ, ϕ соответственно: N, M.

2. Нахождение шагов равномерной сетки по переменным θ, ϕ и узловых

значений переменных θi (i =			), ϕj ( j =		).
		1,N		1,M
3. Вычисление безразмерных параметров: параметра Λ по формуле
Λ = 6μωR2	(δ2 p ) ;
2	o

для схемы А: присвоение параметрам λIi и λIIi (i =1, 5 ) значений относительных

смещений центра ротора и параметров погрешностей геометрии согласно (10.1), (10.2);

для схемы В: нахождение относительных величин εi(k ) (i =1, 3 ; k = I, II) по формулам εi(k ) = ei(k ) δ с учетом (10.10), (10.11).

4. Определение распределения давления слоя жидкости

для схемы А: вычисление коэффициентов сеточных уравнений, полученных на основе дифференциальных уравнений (10.5), (10.6);

применение алгоритма прогонки для решения систем линейных алгебраических уравнений трехдиагональной структуры с целью

определения функций ( X (jk ) )i ,	(Y j(k ) )i (k = I, II; j =		; i =		);
		1, 4		1, n
вычисление узловых значений	функции распределения давления
p(θi, ϕj) ≡ p(( X (jk ) )i , (Y j(k ) )i ,ϕj) по формуле (10.7).

для схемы В: вычисление коэффициентов cs (s =1, m; m = (N – 2)L; L = 2M – 2)

матрицы С и коэффициентов as q (s,q =1, m ) матрицы A матричного

уравнения (10.28) с учетом (10.25), (10.26);

приведение матрицы А к верхнетреугольному виду и нормализация к единице коэффициентов главной диагонали полученной матрицы;

вычисление значений	~(k )	(t =1, m ) методом обратной прогонки;
	pt

нахождение искомых сеточных значений функции распределения

давления

i(kj )

(k = I, II;

i =

;

j =

) путем

p

2, N −1

1, M −1

переобозначения

(k )

~(k )

.

pi j

= pt

174

5. Нахождение масштабного коэффициента для сил	K	F	= πp R2	; значений
			o 2

безразмерных проекций F xi (i =1, 3 ) результирующей реакции подвеса на оси СКxi статора

для схемы А: по формулам (10.8) методом Симпсона численного интегрирования функции одной переменной;

для схемы В: по формулам (10.29) методом Симпсона численного интегрирования функции двух переменных.

6. Вычисление масштабного коэффициента для моментов KМ = πμω R24 δ и

значений безразмерных возмущающих моментов и безразмерного момента сопротивления

для схемы А: по формулам (10.9) с помощью метода интегрирования Симпсона функции одной переменной;

для схемы В: по формулам (10.30) с помощью метода интегрирования Симпсона функции двух переменных.

7.Определение размерных значений проекций результирующих

гидродинамических сил Fxi и моментов MOxi (i =1, 3 ): Fxi = KF F xi и MOxi = KM M Oxi соответственно.

На основе приведенных алгоритмов разработана программа «Анализ погрешностей геометрии сферического гидродинамического подвеса гироскопа (АПГ – СГПГ)» в системе Турбо-Паскаль 7.0, на которую получено Свидетельство об официальной регистрации программы для ЭВМ в РосАПО [26]. Программа предназначена для численного моделирования распределения давления и реакций гидродинамических подвесов, выполненных по схеме А – плоскость смещения усеченных и имеющих полюсные отверстия сферических сегментов, составляющих статор, перпендикулярна оси вращения ротора – и по схеме В – вектор сдвига центров оснований сегментов лежит в плоскости, проходящей через ось вращения ротора.

10.5. Численные оценки влияния погрешностей формы статора на реакции гидродинамического подвеса

миниатюрного шарового гироскопа

Численное моделирование проведено для реальных гидродинамических подвесов с параметрами: давление в камере pо = 5 105 Па; плотность жидкости ρ = 1,75 103 кг/м3; динамический коэффициент вязкости жидкости μ = 7,4 10–4 Па с

при рабочей	температуре	tо = 50 оС;	радиус ротора принимает	значения
R2 = 0,6 10–3 м;	1,185 10–3 м;	1,5 10–3 м;	2,5 10–3 м; номинальное	значение

радиального зазора между статором и ротором при концентричном положении сфер составляет δ = 5 10–6 м; угловая скорость вращения ротора ω = 3,14 103– 18,85 103 рад/с; плотность материала ротора ρр = 8 103 кг/м3.

175

К вопросу левитации ротора в идеальном гидродинамическом подвесе.

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

При условии идеального сферического подвеса (χ = χ3I = χ3II = 0) с параметрами

R2 = 1,185 10–3 м и ω = 6,28 103 рад/с уравновешивание силы тяжести ротора гидродинамической силой достигается при относительном эксцентриситете ε = 5,95 10–5. При перегрузках от переносных ускорений в 10g и 100g величина эксцентриситета ε оказывается соответственно равной 5,95 10–4 и 5,95 10–3.

На рис. 10.9 показана поверхность функции распределения избыточного давления при ε = 5,95 10–4. Поверхность имеет малую несимметричность вдоль

координатной оси θ, которая на приве-

денном графике незаметна. Качественный

вид поверхности не зависит от значений

относительного эксцентриситета, изменяя-

ются лишь максимальные и минимальные

значения избыточного давления:

при ε =5,95 10–5:

max −

=

min −

o

= 9,14 10–5;

p

p

o

p

p

при ε = 5,95 10–3 (перегрузка 100g):

max −

=

min −

o

= 9,14 10–3.

p

p

o

p

p

Отрицательные значения избыточного

o = – 9,14 10–4.

(So):

−

давления появляются при больших отно-

p

p

Рис. 10.9. Поверхность

сительных

эксцентриситетах.

Так,

функции распределения

например, ε > 0,49 для рассматриваемого

избыточного давления

миниатюрного ротора при перегрузках,

превышающих 1000g.

Безразмерный коэффициент экваториальной жесткости для подвеса с указанными параметрами имеет значение 2,072. Масштабный коэффициент для перехода к размерному значению коэффициента экваториальной жесткости есть

Kg = 4,412 105 Н м.

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

Принимается, что центры OxI и OxII сферических сегментов, составляющих

каркас статора не совпадают:		OxI OxII		= χ (см. рис.	10.2а, б). За полюс Ox СКxi
статора выбран центр	OxI сегмента			I. Центром		подвеса	назовем точку O*,
занимающую среднее	положение между точками					OxI и	OxII : OxO* = 0,5 χ

(рис. 10.10а, б). Тогда положение центра O ротора в СКxi статора определяется векторным равенством

176

e = 0,5 χ + e* ,

где e* – вектор смещения центра O ротора относительно центра подвеса O*.

		χ	0,5		cos
x* x 1	φ	* 2 Ox , OxI
						I
1		O		O	Ox , Ox
	Ox	e	Ox		(э) O	e (э)
		χ e* O		e			e (э) = ex x	+ ex x
		m р g	cos
					m р g		e(э) = ex x + ex x
		x	x
		x 2					x
		x *2				x

а) для схемы А

б) для схемы В

Рис. 10.10. Векторы смещений в плоскости Оx x1x2

Проекции относительного эксцентриситета ε в рассматриваемой задаче записываются соответственно для схем А и В:

ε	= 0,5χ cosφ+ ε*	,	ε	2	= 0,5χ sin φ+ ε*				,	ε	3	= ε*	;
1	1							2				3
ε	= 0,5χ cos φ+ ε*	,	ε	2	= ε*	, ε	3	= 0,5χ sin φ+ ε* .
1	1				2							3

В случае схемы А зависимости проекции Fx1 гидродинамической силы от проекции ε2 относительного эксцентриситета показаны на рис. 10.11.

Рис. 10. 11. Зависимости гидродинамической реакции Fx1 от проекции ε2 относительного эксцентриситета

Графики построены при различных значениях радиуса ротора R2, угловой скорости ω и заданных величинах сдвига χ центров сегментов, составляющих

177

каркас статора. Зависимости Fx2 (ε1 ) при тех же параметрах подвеса аналогичны.

Величины проекций Fx1 , Fx2 результирующей гидродинамической силы

возрастают при увеличении радиуса ротора и угловой скорости его вращения.

Значения проекций Fx1 , Fx2

не зависят от величины сдвига χ центров сегментов

и угла ориентации φ вектора χ.

Диаграмма расположения векторов результирующих гидродинамических

реакций подвеса F и соответствующих векторов смещения e* центра O ротора

относительно центра подвеса O* изображена на рис. 10.12. На диаграмме векторы

F

x 2*

сил

рассматриваются

приложенными

к

2

F1

центру О ротора

для

всех представленных

положений, также показаны соответствую-

F3

e*3

F5

щие векторы эксцентриситета e

относитель-

e*2

но

полюса

Ox

статора. Для

подвесов

с

O*

e*1

x 1*

параметрами:

R2 = 1,185 10–3 м;

e2

e5*

x 2

e3

e4*

ω = 15,7 103 рад/с; ε* = 0,001 и R2 = 1,5 10–3 м;

e4

e5

F4

ω = 9,42 103 рад/с;

ε* = 0,0001

модули

φ

Ox

χ

x 1 e1

гидродинамических сил F равны 1,15 10–2 Н

и 1,77 10–3 Н соответственно.

2

С помощью такой диаграммы можно оп-

Рис. 10.12. Векторная диаграмма

ределять равновесное положение центра ро-

гидродинамических сил F

тора с учетом перегрузки от переносных

ускорений.

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

смещений ε* в задаче левитации без перегрузки

( mрg = F ) для схемы А зависят от радиуса ротора R2 и угловой скорости ω, как

приведено в табл. 10.1.

Таблица 10.1

R2, м	1,185·10–3	1,5·10–3	2,5·10–3
ω, рад/c
3,14·103	–	–	ε*= 1,13·10–4
6,28·103	ε*= 1,19·10–4	ε*= 9,41·10–5	ε*= 5,65·10–5
9,42·103	ε*= 7,94·10–5	ε*= 6,28·10–5	–
1,26·104	ε*= 5,95·10–5	ε*= 4,70·10–5	–
1,57·104	ε*= 4,76·10–5	ε*= 3,41·10–5	–

Данные в таблице дают представление о величинах относительных смещений ε*, которые оказываются весьма малыми. При перегрузках от переносных ускорений в 10g и 100g величины ε* возрастают в 10 и 100 раз соответственно.

178

В подвесе, представленном схемой В, вектор F результирующей гидродинамической силы имеет не только экваториальную F (экв) = Fx1 x1 + Fx2 x2 ,

но и осевую F x3 составляющие – если величина сдвига χ ≠ 0. Зависимости экваториальных Fx1 , Fx2 и осевых Fx3 проекций гидродинамических сил от

проекций ε*i (i =1, 3) относительного эксцентриситета ε* показаны на рис. 10.13,

10.14 для подвеса с параметрами R2 = 1,185 10–3 м, ω = 15,71 103 рад/с и указанными величинами сдвига χ центров сегментов.

Графики 1, 2: Fxi (ε1* ) ; 3, 4: Fxi (ε*2 ) ; 5, 6: Fxi (ε*3 ) ; i = 1, 2.

Рис. 10.13. Зависимости экваториальных проекций Fx1 , Fx 2 гидродинамических сил от проекций ε*i (i = 1, 3) относительного эксцентриситета

Графики 1, 4: Fx3 (ε1* ) ; 2, 4: Fx3 (ε*2 ) ; 3, 4: Fx3 (ε*3 ) .

Рис. 10.14. Зависимости осевых проекций Fx3 гидро-

динамических сил от проекций ε*i (i = 1, 3) относительного эксцентриситета

179

Если

центры сегментов, составляющих каркас

статора, не совпадают

(χ ≠ 0),

то

значения

проекций

Fxi

(i =

)

главного

вектора

F

1, 3

гидродинамических сил

зависят

от

экваториальных ε* ,

ε*

и осевой

ε*

составляющих вектора смещения ε* .

1

2

3

Величины осевых проекций

вектора

F

значительно

меньше

значений

экваториальных

проекций:

F

≈10−3 F

x1

,

−3 F

x3

F

≈10

. Осевая составляющая

F

является практически постоянной, и

x3

x2

x3

осевая жесткость не обнаруживается.

На рис. 10.15а, б изображены диаграммы взаимного расположения экваториальных составляющих гидродинамических сил F (экв) и экваториальных

смещений e* (э) центра O ротора относительно центра подвеса O*, а также составляющих e (э) эксцентриситетов e .

Ox

3(э)

1(э)

2(э)

e

e

e

2*(э)

4(э)

e

cos

0,5

χ

e

2(экв)

(э)

F

*

e3*

O*

(экв)

x 2, x 2

e1*(э)

F3

(э)

F (экв)

e

4*

1

F

(экв)

4

x 1, x *1

e i(э) = e *i (э) (i = 1, 4)

(экв)

F

2

2(э)

e

x 2, x 2*

O*, Ox

e 3(э)

(э)

e1

(экв)

4(э)

F1

e

F

(экв)

3

F

(экв)

4

x 1, x *1

а) χ = 10–7 м; φ= π/4 б) χ = 10–7 м; φ = π/2

Рис. 10.15. Векторные диаграммы экваториальных составляющих гидродинамических сил F (экв)

Векторные диаграммы (см. рис. 10.15а, б) соответствуют гидродинамическому

подвесу с	параметрами R2 = 1,185 10–3 м; ω = 15,71 103 рад/с; ε* = ε* (э) = 10–4;
χ = 10–7 м. При ориентации φ = π/4 сдвига центров сегментов (см. рис. 10.15а)
направления	экваториальных составляющих	F (экв) главных векторов

гидродинамических сил почти совпадают во всех положениях центра О ротора

относительно		центра О* подвеса. Величины экваториальных составляющих
F (экв) = F 2	+ F 2		одинаковы и равны ≈ 0,164 Н. При этом значения осевых
x1		x2

3

проекций Fx3 ≈ 1,68·10–4 Н. Отклонения величин F = ∑Fxi2 , вычисленных при

i=1

различных направлениях вектора e* (э), от их среднего значения составляют не более 0,67%. Полученные гидродинамические силы F уравновешивают силу инерции ротора от переносных ускорений ≈ 298g.

180