Введение. О проблеме транспортных перевозок в северных районах россии

В освоении арктических регионов РФ значительна роль транспортной проблемы [1]. Если в среднем по стране доля транспортных расходов в себестоимости продукции составляет 10…12%, то на Севере она может достигать 30…60% и более. Это связано с природно-климатическими и экономическими особенностями, определяющими особый подход к формированию транспортной системы:

• низкая плотность населения (до 0,3 чел./км²) и большие расстояния между малыми населенными пунктами, что определяет небольшие пассажиропотоки.
• очаговая деятельность при освоении сырьевых месторождений;
• суровые климатические условия (годовой ход температур – до 100°), короткое лето, полярная ночь;
• вечная мерзлота, значительное количество болот, рек, озер, экологическая уязвимость экосистем;
• повышенные цены на топливо, оборудование и товары из-за разового завоза и хранения межсезонных запасов;
• высокая стоимость создания и эксплуатации транспортной инфраструктуры, слабое развитие наземных транспортных коммуникаций.

Местные авиационные перевозки дорогостоящи; морское и речное судоходство возможны в короткий летний период; применение автотранспорта затрудняется, в т.ч., по экологическим причинам, из-за сезонности, и др.

Прорыв может быть достигнут применением судов на воздушной подушке (СВП) в составе существующей транспортной системы, так как они практически безальтернативны для круглогодичных перевозок как вне, так и внутри Арктики [1]. Амфибийность позволяет им выходить на пологий берег и преодолевать большие расстояния по слабо пересеченной местности. При частых (каждые 20…50 км) остановках экономические показатели СВП могут существенно превосходить показатели самолетов и вертолетов.

Однако, технические решения, лежащие в основе конструкций распространенных в настоящее время СВП, не отвечают требованиям, предъявляемым указанными условиями эксплуатации. Эффективность СВП резко снижается на пересеченной местности, при движении по торосам, на уклонах и косогорах. Они неустойчивы и не имеют необходимой для движения по суше, извилистым узким рекам и в иных случаях стесненной обстановки траекторной управляемости. Наконец, даже на небольших пологих уклонах (более 5^о) такие СВП плохо удерживаются на траектории. Потеря подвижности СВП на местности, как правило, происходит из-за потерь воздуха через зазоры между гибким ограждением (ГО) и подстилающей поверхностью. В некоторых СВП этот эффект отчасти удается компенсировать за счет мероприятий, увеличивающих сопротивление движению из-за трения ГО о грунт. С этим, однако, связан повышенный износ ГО и риск его повреждения, а также рост сопротивления движению. Кроме того, классическая конструкция ГО содержит внутренние перегородки, доступ к которым при выключенных нагнетателях и осевшем СВП практически отсутствует, что делает невозможным их ремонт в полевых условиях. Это может вызывать большие трудности при эксплуатации СВП в малонаселенных районах.

Анализ условий эксплуатации позволил сформировать основные требования к наиболее массовым транспортным СВП для указанных условий эксплуатации:

1. Грузоподъемность не менее 10 т (приспособленность к размещению стандартного контейнера ISO, а также возможность сквозной погрузки на грузовую платформу);
2. Способность двигаться по всем видам подстилающих поверхностей со скоростью не менее 40…45 км/ч;
3. Габаритные размеры, обеспечивающие беспрепятственное движение по рекам низших технических категорий;
4. Способность уверенно передвигаться над неровными подстилающими поверхностями (с высотой или глубиной неровностей до 1,6…1,8 м, близких по длине к длине аппарата) или над водными поверхностями (высота волны до 1,8 м);
5. Способность создавать уровень воздействия на подстилающую поверхность, который исключает ее невосстановимые повреждения;
6. Способность преодолевать затяжные подъемы крутизной до 15^о., а также устойчиво двигаться по склонам крутизной до 10^о;
7. Способность к сохранению управляемости и курсовой устойчивости при ветре со скоростью до 10…15 м/с под всеми курсовыми углами.

В [2–9] показано, что этот комплекс требований может быть реализован в СВП, имеющих контактные движители (колеса) с частичной разгрузкой при помощи воздушной подушки. Основная масса поддерживается ВП, а нагрузка на контактный движитель регулируется таким образом, чтобы обеспечить уровень тягово-сцепных свойств, необходимый для уверенного маневрирования, преодоления подъемов и удержания машины на курсе при боковом воздействии.

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

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

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

Комплексная математическая модель СВП

Гибкое ограждение классического типа

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

Истечение через щель считается пропорциональным площади зазора под элементом и квадратному корню удвоенного избыточного давления в подушке, отнесенного к плотности воздуха. Так же вводится эмпирический поправочный коэффициент, рассчитываемый в общем случае как \(\mu_i=E(\theta,\tau/h,\tau/b)\mu_0(\theta)\), где \(E\) – коэффициент, связанный с податливостью подстилающей поверхности и геометрией навесных элементов, \(\mu_0\) – коэффициент расхода воздуха для ГО в виде сплошной полосы над твердым экраном, \(\theta\) – угол наклона навесных элементов или ГО к горизонту, \(\tau\) – ширина щели под ГО, \(b\) – ширина навесных элементов, \(h\) – глубина впадины под ГО, образованной податливой поверхностью. Для жесткого экрана \(E(\theta,\tau/h,\tau/b) = 1\).

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

В объеме под ГО изменение давления рассчитывается решением следующего дифференциального уравнения для каждой секции [11]:

\[\frac{\partial P}{\partial t}=\frac{c^{2}}{V}(G_{c}- \sum_{i=0} ^{N-1} Q_{i}\rho_{a})\]

Здесь \(c\) – скорость звука в воздухе, \(P\) – давление в секции подушки, \(V\) – объем секции, \(G_c\) – массовый расход на выходе нагнетателя, \(Q_i\) – объемный расход в \(i\)-том элементе (всего \(N\) элементов), \(\rho_a\) – плотность воздуха. Для текущего приближения считается, что подушка снабжена четырьмя нагнетателями, по одному на каждую секцию.

Расход нагнетателя связан с давлениями на выходе из нагнетателя и в объеме подушки и геометрией объема подушки [11] следующим дифференциальным уравнением:

\[\frac{\partial G_{c}}{\partial t}=\frac{A}{L}\left(p_{c}-P\right)\]

Здесь \(G_c\) – массовый расход на выходе нагнетателя, \(A\) – характерная ширина потока воздуха, \(L\) – характерная длина потока воздуха, \(p_c\) – давление на выходе из нагнетателя, \(P\) – давление в секции подушки.

Сопротивление движению

Сопротивление движению судна на воздушной подушке определяется как сумма следующих аэродинамических и гидродинамических эффектов [12].

Импульсное сопротивление связано с ускорением поступающего в нагнетатель воздуха до скорости движения судна. Оно пропорционально массовому расходу через нагнетатели и скорости набегающего на судно воздушного потока.

Профильное сопротивление – это сопротивление, обусловленное аэродинамическим сопротивлением корпуса, и рассчитывается по стандартным формулам для турбулентного обтекания исходя из формы корпуса [13].

Сопротивление струй вызвано истечением воздуха из-под воздушной подушки с учетом закона сохранения импульса:

\[\mathbf{R_{str}}=\rho_{a}\sum_{i=0}^{N-1}\frac{Q_{i}^{2}}{\tau_{i}L_{i}}\mathbf{n_i}\]

Здесь \(\rho_a\) – плотность воздуха, \(Q_i\) – объемный расход воздуха из \(i\)-го элемента, \(\tau_i\) – ширина зазора под элементом, \(L_i\) – длина нижней кромки данного элемента, \(\mathbf{n_i}\) – нормаль к огибающей ГО в середине элемента.

Волновое сопротивление связано с работой на образование впадины в воде в пространстве под воздушной подушкой:

\[\begin{aligned} R_{wave} & =\kappa_{w}^{1}\frac{G_{L}Mg}{\alpha_{\Pi}}\\ \alpha_{\Pi} & =\frac{S}{LB}\\ G_{L} & =\frac{P_{avg}}{\rho_{w}gL}\end{aligned}\]

Здесь \(\alpha_\pi\) – коэффициент заполнения площади подушки, \(G_L\) – продольный коэффициент загрузки, \(P_{avg}\) – среднее давление в подушке, \(\rho_w\) – плотность воды, \(g\) – ускорение свободного падения, \(L\) – длина подушки, \(B\) – ширина подушки, \(S\) – площадь подушки, \(M\) – масса судна, \(\kappa_w^1\) – эмпирический коэффициент, задаваемый графически [14].

Другие составляющие гидродинамического сопротивления (от брызг, касания ГО воды и др.) принято обозначать как остаточное сопротивление \(R_{res}\) [15]:

\[R_{res}=\frac{c_{\Gamma0}L}{\left(\tau/L\right)^{0,34}}\sqrt{S}\frac{\rho_{w}v^{2}}{2}+\left[\frac{5,09}{\left(G_{L}\rho_{w}g\right)^{0,259}}-1\right]R_{wave}.\]

где \(\rho_w\) – плотность воды, \(v\) – скорость судна, \(\tau\) – средняя величина зазора, \(L\) – длина подушки, \(S\) – площадь подушки, \(c_{ГО}=15.4 \cdot 10^{-6}\), \(G_L\) – продольный коэффициент загрузки, \(g\) – ускорение свободного падения, \(R_{wave}\) – волновое сопротивление.

Нагнетатели

В основе описания рабочих процессов в аксиальном нагнетателе лежит полуэмпирический подход, описанный в [16–20].

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

Поэтому при разработке данной модификации математической модели необходимо было как можно точнее из основных принципов описывать рабочие процессы, не выходя за рамки неплохо обоснованного канального приближения. Применено теоретическое описание газодинамических процессов на основе балансовых соотношений [21–23] , отражающих сохранение массы, импульса и энергии. Требования описания многих режимов не позволяют использовать упрощающие предположения, с успехом применяющиеся при проектировании [17; 18].

Здесь и далее скорости в лабораторных координатах обозначены \(c\), скорости во вращающихся координатах, связанных с рабочим колесом – \(w\), скорость вращения колеса в рассматриваемом месте колеса – \(u\), скорость звука – \(a\) или \(v_s\), индексом \(u\) или \(\theta\) обозначены проекции скоростей на направление азимутального вращательного движения в данном месте колеса, индексом \(a\) – проекции скоростей на направление оси симметрии, индексом \(r\) – радиальные составляющие скоростей. Индексом \(1\) отмечены параметры на входе в подвижные лопатки, \(2\) – на выходе из подвижных рабочих лопаток, \(3\) – в выходном патрубке нагнетателя (после успокоения потока).

Рассмотрим описание физической и математической постановки для модуля, предназначенного для расчета процессов в нагнетателе.

Полагалось, что на входе в нагнетатель известны давление \(p_1\), температура газа \(T_1\) и массовый расход \(G\). Скорость потока на входе параллельна оси вращения колеса. По известным давлению и температуре на входе в подвижные лопатки (состояние 1) находятся значения остальных параметров: плотности, скорости потока, скорости звука.

\[\rho_1 = \frac{m_{mol} p_1}{k_B T_1}\]

\[c_{1a} = \frac{G}{\rho_1 F_0}\]

\[c_1 = c_{1a}\]

\[c_{1\theta} = 0\]

\[v_{s1} = \sqrt{\frac{k k_B T_1}{m_{mol}}}\]

\[M_1 = \frac{c_1}{v_s1}\]

где \(F_0\) – площадь сечения каналов на входе в лопатки, перпендикулярного оси вращения.

По найденным характеристикам потока, набегающего на рабочее колесо, определяются модуль и направление относительной скорости \(w_1\), число Маха набегающего потока во вращающихся координатах рабочего колеса \(M_{w1}\), угол входа потока в каналы рабочего колеса \(\beta_{11}\), температура торможения \(T_{w1}^*\):

\[u = \omega r_0\]

\[w_{1\theta} = c_{1\theta} - u\]

\[w_{1a} = c_{1a}\]

\[\beta_{11} = \arctan\left(\frac{w_{1a}}{w_{1\theta}}\right)\]

\[w_1 = \sqrt{w_{1a}^2 + w_{1\theta}^2} M_{w1} = \frac{w_{1t}}{v_{s1}}\]

\[T_{w1}^* = T_1 + \frac{w_1^2}{2C_p}\]

где \(\omega\) – угловая скорость, \(C_p\) – удельная теплоемкость газа.

Находится эффективное критическое сечение канала и по нему – число Маха на выходе и остальные параметры в сечении “2” как в координатах колеса, так и в лабораторных координатах. При этом возникает необходимость решения трансцендентного уравнения:

\[F_{21} = F_0 \sin \beta_{11}\]

\[F_{2cr} = F_{21} \frac{ M_{w1} \left( \frac{1}{2} k + \frac{1}{2} \right) ^{\frac{1}{2}\frac{k+1}{k-1}} }{ \left(1+\frac{1}{2} M_{w1}^2(k-1)\right) ^{\frac{1}{2}\frac{k+1}{k-1}} }\]

Оно решалось численно, результаты аппроксимировались и применялись при дальнейших расчетах.

\[M_{w2} = M'\frac{F_2}{F_{2cr}}\]

\[c_{cr} = \sqrt{\frac{w_1^2 (k-1)+2 v_{s1}^2}{k+1}}\]

\[T_{cr} = \frac{m_{mol} c_{cr}^2}{k k_B}\]

\[T_2 = T_{cr} \frac{1+\frac{k-1}{2}}{1+\frac{k-1}{2} M_{w2}^2}\]

\[a_2 = \frac{k k_B T_2}{m_{mol}}\]

\[w_2 = M_{w2} a_2\]

\[w_{2a} = w_2 \sin \alpha_2\]

\[w_{2\theta} = w_2 \cos \alpha_2\]

\[\rho_2 = \frac{G}{w_{2a} F_0}\]

\[p_2 = \frac{\rho_2}{m_{mol} k_B T_2}\]

\[c_{2a} = w_{2a}\]

\[c_{2\theta} = w_{2\theta} − u\]

\[c_2 = \sqrt{c_{2a}^2 + c_{2\theta}^2}\]

\[M_2 = \frac{c_2}{a_2}\]

где \(k_B\) – постоянная Больцмана, \(m_{mol}\) – масса молекулы газа, индекс \(cr\) относится к параметрам в эффективном критическом сечении.

По найденным параметрам в сечении “2” находится давление торможения в предположении, что восстанавливается осевая составляющая скорости, а азимутальная диссипирует:

\[p_{3i} = p_2 + \frac{ρ_2 c_{2a}^2}{2}\]

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

\[ Y = \frac{p^*_{in} - p^*_{out}}{p^*_{out} - p_{out}}\]

где \(p^*_{in}\), \(p^*_{out}\) – полное давление на входе и выходе, \(p_{out}\) – статическое давление на выходе лопаток.

Значение \(Y_P\) находится из соотношения

\[Y_P = \mu \chi \left(Y_{Pr} + \left(\frac{\alpha_{bin}}{\alpha_{bout}}\right)^2 (Y_{Pi} - Y_{Pr})\right)\cdot\left(\frac{5t}{C}\right)^\frac{\alpha_{bin}}{\alpha_{bout}}+Y_S\]

где значения \(Y_{Pr}\) и \(Y_{Pi}\) снимаются с графиков [24] в зависимости от угла лопатки на выходе \(\alpha_{bout}\) и от отношения шага между лопатками \(s\) к длине лопатки \(C\); здесь \(\alpha_{bin}\) – угол лопатки на входе потока, \(t\) – толщина лопатки, обычно составляющая

\[t=(0,15..0,25) C\]

Влияние вязкости учитывается фактором \(\chi\), который зависит от числа Рейнольдса \(Re = w_2 \frac{C}{\nu}\):

при \(Re < 2 \cdot 10^5\), \(\chi = \left(\frac{Re}{2\cdot10^5}\right)^{-0,4}\)

при \(2\cdot 10^5 \leq Re \leq 10^6\), \(\chi = 1\)

при \(Re > 10^6\), \(\chi = \left(\frac{Re}{10^6}\right)^{−0,2}\)

Фактор \(\mu\) применяется для учета вторичных вихрей

\[\mu = 1 - \frac{Z_{TE}}{h}\]

\[Z_{TE} = \left(\frac{0,1 F_\theta^{0,79}}{\sqrt{C_R}\left(\frac{h}{C}\right)^{0,55}}+32,7 \frac{\delta^*}{h}\right)h\]

Здесь \(h\) – ширина лопатки в радиальном направлении, \(\delta^*\) – толщина погранслоя, она составляет \((0,008..0,06) h\), фактор

\[F_\theta = \frac{2s}{C_x}(\tan \alpha_{in}+\tan \alpha_{bout})\cos^2 \alpha_m\]

\[\alpha_m = \frac{\alpha_{bin} + \alpha_{bout}}{2}\]

здесь \(\alpha_{in}\) – угол входа потока, \(C_x\) – длина проекции лопатки на ось.

Учет вклада вторичных потерь:

при \(0\leq x''\leq 0,3\)

\[Y_S = Y_{S0}\left(\exp(0,9 x'')+13x''^2+400x''^4\right)\] Y_S =Y_S0(exp(0,9x”)+13x”²+400x”⁴),

при \(-0,4 \leq x'' \leq 0\)

\[Y_S =Y_{S0} \exp(0,9 x'')\]

\[x'' = \frac{\alpha_{in}-\alpha_{bin}} {\pi-(\alpha_{bin}+\alpha_{bout})} \left(\frac{\cos \alpha_{bin}}{\cos \alpha_{bout}}\right)^{-1,5} \left(\frac{d_{in}}{C}\right)^{-0,3} \]

Учет потерь из-за нерасчетного угла входа проводится с помощью полуэмпирического соотношения [25] для связанного с этим изменения \(\Delta \phi_P^2\) квадрата отношения \(\phi = {c_2}/{c_{2\tau}}\) (\(c_{2\tau}\) – скорость без учета этого вида потерь):

при \(x \geq 0\)

\[\Delta \phi_P^2 = 3,711\cdot 10^{-7} x^8 - 5,318\cdot 10^{-6} x^7 + 1,106\cdot 10^{-5} x^6 + 9,017\cdot 10^{-5} x^5 - 1,542\cdot 10^{-4} x^4 - 2,506\cdot 10^{-4} x^3 + 1,327\cdot 10^{-3} x^2 - 6,149\cdot 10^{-5} x \]

при \(x < 0\)

\[\Delta \phi_P^2 = 1,358 \cdot 10^{-4} x - 8,720\cdot10^{-4}x \]

\[ x = \left(\frac{d_{in}}{s}\right)^{−0,05} \left(\arccos\frac{C_x}{C}+\alpha_{bin}\right)^{−0,2} \left(\frac{\cos \alpha_{bin}}{\cos \alpha_{bout}}\right)^{−1,4} (\alpha_{in}−\alpha^*_{in}) \]

здесь \(\alpha^*_{in} \approx a_{bin}\) – расчетный угол входа потока, \(d_{in}\) – диаметр закругления лопатки на входе.

Описанная комплексная математическая модель должна модифицироваться в ходе разработки при уточнении и принятии новых технических решений; так, может производиться пересчет в подробных 3D постановках с учетом изменения параметров конструкции, режимов и др.

Винтовые движители

Методика расчета фактически разработана на основе классической теории винта [???; 26–34]. Эта теория часто дает более надежные данные, чем математическое моделирование (в т.ч. в трехмерных постановках на основе методов низкого порядка точности при числе расчетных ячеек порядка миллиона), поэтому она не теряет актуальность до сих пор [35].

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

В данном случае, в отличие от классического подхода, циркуляция и скорость набегающего на лопатки потока заранее не известны. Это существенно меняет алгоритм реализации. Тангенциальная \(u_1(\Gamma) = \frac{k\Gamma}{4\pi r}\) и осевая \(v_1(\Gamma) = -V/2 + \sqrt{V^2/4 + u_1 (\omega r - u_1) + 2 \int_r^R u_1^2/r^3 dr}\) составляющие индуктивной скорости в плоскости винта представляют собой зависимости от циркуляции. Затем находится скорость набегающего потока \(W_1(\Gamma) = \sqrt{(\omega r - u_1(\Gamma))^2+(V+v_1(\Gamma))^2}\), рассчитывается циркуляция \(\Gamma = 0,5 b C_y W_1(\Gamma)\). Полученное уравнение решается методом бисекции относительно \(\Gamma\). Расчет начинается с концевых сечений винта. Здесь \(k\) – число лопастей, \(r\) – текущий радиус, \(\Gamma\) – циркуляция на данном радиусе, \(V\) – скорость осевой обдувки винта, \(\omega\) – угловая скорость вращения винта, \(R\) – наибольшее значение радиуса винта, \(b\) – длина хорды лопасти. Интеграл под корнем в выражении для \(v_1(\Gamma)\) соответствует центробежным силам, которыми в данном случае можно пренебречь, благодаря наличию кольцевой насадки вокруг винта, т.к. она ограничивает движение в радиальном направлении. По этой же причине не учитываются концевые потери: вихревые линии не уходят в область свободного потока, а кончаются на кольцевой насадке.

За винтом осевая и тангенциальная составляющие скорости в два раза больше, чем в плоскости вихревого диска:

\[u_2=2 u_1\]

\[v_2=2 v_1\]

Для области \(r>R\) индуктивная скорость полагается равной нулю.

Аэродинамические коэффициенты \(C_y\), \(C_x\) определяются профилем винта. Они известны для ряда профилей из подробных экспериментальных данных [36, p. [@lushin1994atlas]]. Эти величины находятся с помощью линейной интерполяции и экстраполяции по табличным данным для профиля. В данном примере взят профиль СибНИА С-11.

Если полагать плотность воздуха постоянной в окрестностях винта и равной \(\rho\), подъемная сила, действующая на элемент лопасти и направленная перпендикулярно скорости набегающего потока, может быть найдена, как \(dR_1 = \rho\Gamma W_1 dr\), а сила профильного сопротивления, направленная вдоль скорости набегающего потока, как \(dR_p=\frac{C_x}{C_y}dR_1\).

Сила тяги \(dT_1\) и сила сопротивления \(dQ_1\) для элемента лопасти могут быть найдены, как проекции суммарной силы \(dR_1+dR_p\) на ось вращения и плоскость вращения соответственно.

Мощность \(dL\), необходимая для вращения элемента лопасти, лежащего на окружности радиуса r, может быть определена из момента силы сопротивления \(dM_{Q1}= −dQ_1 r\):

\[dL=dM_{Q1} \omega\]

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

Нередко винты АСВП располагаются частично в зоне аэродинамической тени от корпуса. Основное воздействие от этого на рабочие процессы связано с тем, что скорость воздуха на входе \(V\) не постоянна по площади винта, а зависит от вертикальной и поперечной координат. Распределение скорости \(V(x,y)\) в первом приближении можно рассчитать отдельно по обтеканию корпуса без учета влияния винтов. Тогда взаимодействие лопасти винта и набегающего потока в данный момент зависит от углового положения лопасти. Для незатененных участков расчет полностью аналогичен приведенному выше, а для затененных – проводится с переменной локальной скоростью \(V(x,y)\), при этом получаются различные значения углов потока относительно лопаток, сил, моментов, выходных скоростей и пр. в зависимости от углового положения лопасти.

В качестве начального приближения можно задать \(V(x,y)\) в виде ступенчатой функции. Тогда можно просчитать по представленной выше методике обтекание лопатки при скоростях, соответствующих значениям скорости на уровнях “ступенек”. При задании аэродинамической тени с помощью \(N\) уровней скорости расчет проводится соответственно \(N\) раз. Эти результаты используются в соответствующих участках заметаемого винтом круга при интегрировании по этому кругу сил, моментов, определении выходных скоростей и пр. В качестве начального варианта выбран двухуровневый вариант.

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

Рулевой комплекс

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

Для расчета сил и моментов используется система координат: ось \(OZ\) направлена по оси винта, ось \(OY\) – параллельно прямым образующим рулей, ось \(OX\) – так, чтобы система координат была правой. Центр \(O\) находится на оси вращения винта.

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

Поперечная сила, действующая на элемент руля, аналогичная подъемной силе крыла. Сила равна

\[dR_1 = \frac{1}{2} C_v \rho v^2 b_v dy\]

где \(C_v\) – аэродинамический коэффициент подъемной силы для профиля руля, \(b_v\) – длина хорды руля, \(v\) – локальная скорость набегающего на руль потока.

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

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

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

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

Аэродинамические коэффициенты в принципе определяются так же, как и для винта.

Моменты сил определяются по векторной формуле

\[d \mathbf{M}_{R1} = \mathbf{r}\times d\mathbf{R}_1\] \[d\mathbf{M}_{Rp} = \mathbf{r} \times d\mathbf{R}_p\]

Для получения суммарных сил и моментов, действующих на решетку рулей, производится интегрирование по всем рулям, реализованное в виде суммирования по элементам рулей.

Описанная методика позволяет рассчитать в базовом приближении все основные характеристики движителя и рулей − тягу и управляющий момент, приложенные к АСВП, момент сопротивления на валу движителя − в зависимости от управляющих воздействий (поворот лопаток движителя, поворот рулей, частоты вращения вала движителя) и внешних условий (скорость ветра, скорость АСВП, плотность воздуха).

Методика расчета конического элемента гибкого ограждения

Нам не известны работы по теории конического элемента гибкого ограждения, по которым можно было бы построить математическую модель такого типа ГО. Это заставило провести собственное численное моделирование. Более того, анализ известных из литературы методов расчета гибких оболочек показал их непригодность для применения в составе комплексной математической модели – закрытость кодов, невозможность в одном методе учесть сопротивление изгибу и сдвигу, конечное растяжение, возможные жесткие элементы, растяжки, сложность задания геометрии ГО, большое время счета и др. Поэтому нами применялась следующая методика.

Нерастяжимую мнущуюся поверхность представим в виде регулярной сетки, близкой к прямоугольной. В каждом ее узле сосредоточена определенная масса (в общем случае, массы узлов могут различаться). Узлы соединены между собой невесомыми пружинами (ребрами). Искомая коническая поверхность имеет \(N_z+1\) слоев, каждый из которых состоит из \(N_\phi\) узлов. Положение и скорость \((i, j)\)-ого узла определяется векторами \(\mathbf{x}_{i,j}\) и \(\mathbf{v}_{i,j}\).

Динамика каждого узла оболочки во времени определяется уравнениями движения Ньютона. В данной модели на каждый узел, представляемый материальной точкой, действуют силы тяжести, растяжения/сжатия, сопротивления изгибу, трения, давления и реакции опоры.

Сила тяжести \(\mathbf{F}_g = \rho S h \mathbf{g}\), где \(\rho\) – плотность материала, \(S\) – площадь ячейки, \(h\) – толщина оболочки, \(\mathbf{g}\) – вектор ускорения свободного падения.

Сила растяжения/сжатия \(\mathbf{F}_t = \sum_j k_t (1-\frac{L_j}{d_j})\mathbf{d}_j\), где суммирование ведется по всем смежным ячейкам (в том числе по диагональным направлениям), \(k_t\) – коэффициент упругости (для диагональных направлений коэффициент упругости считается равным \(0,3 k_t\)), \(\mathbf{d}_j\) – вектор, соединяющий текущую ячейку с \(j\)-той ячейкой, \(L_j\) – равновесная длина связи.

Сила реакции опоры \(\mathbf{F}_{nr} = k_{nr} l^2 \mathbf{n}\), где \(k_{nr}\) – коэффициент жесткости взаимодействия, \(l\) – глубина проминания подстилающей поверхности, \(\mathbf{n}\) – единичный вектор нормали к поверхности.

Сила трения о воздух \(\mathbf{F}_{fr}^{air} = - k_{fr}^{air} \mathbf{v}\), где \(k_{fr}^{air}\) – коэффициент торможения о воздух, \(\mathbf{v}\) – скорость узла относительно воздуха.

Сила трения о подстилающую поверхность \(F_{fr} = - k_{fr} F_{reaction} \mathbf{\tau},\) где \(k_{fr}\) – коэффициент трения о землю, \(\mathbf{v}\) – скорость ячейки относительно подстилающей поверхности, \(\mathbf{\tau} = \frac{\mathbf{v} - (\mathbf{v} \cdot \mathbf{n})\mathbf{n}}{\left|\mathbf{v} - (\mathbf{v} \cdot \mathbf{n})\mathbf{n}\right|}\) – единичный вектор в направлении проекции скорости на подстилающую поверхность.

Сила давления \(\mathbf{F} = P S \mathbf{n}\), где \(\mathbf{n}\) – единичный вектор внешней нормали к поверхности оболочки, \(P\) – избыточное давление под оболочкой, \(S\) – площадь ячейки. Сила сопротивления изгибу на каждом узле вычисляется для вертикального и горизонтального направлений. Изгибающий момент вычисляется по формуле

\[M = (C - C_0) E J\]

где \(E\) − модуль Юнга, например, для резины \(E \approx 10^7 \text{Па}\); \(J\) − момент инерции сечения в соответствующем направлении. Для плоского сечения \(J = \frac{L h^3}{12}\), где \(h\) − толщина гибкой оболочки, \(L\) − средняя длина двух противолежащих ребер (в горизонтальном или вертикальном направлении), исходящих из данного узла, \(C\) – задаваемая координатами узлов шаблона кривизна поверхности в данной точке (в горизонтальном или вертикальном направлении), \(C_0\) – равновесные значения кривизны, при которых силы не возникают.

В случае близкой к равномерной сетки, \(|C|= \frac{1}{R}\), \(R = \frac{a}{2\sin\alpha}\), \(a = 2 L \sin\frac{\alpha}{2}\), \(\cos\alpha = \frac{\mathbf{L}_1 \cdot \mathbf{L}_2}{L_1 L_2}\), \(L = \frac{L_1+L_2}{2}\), где \(R\) – радиус кривизны поверхности, \(\alpha\) – угол между двумя противолежащими ребрами, \(\mathbf{L}_{1,2}\) – векторы противолежащих ребер.

Силы рассчитываются для трех узлов следующим образом:

\(\mathbf{F}_{1,2} = -\frac{M}{L_{1,2}} \mathbf{n}_{1,2}\) – для узлов в конце \(\mathbf{L}_{1,2}\) соответственно,

\(\mathbf{F} = -(\mathbf{F_1}+\mathbf{F_2})\) – для среднего узла.

Здесь \(n_{1,2} = \frac{\mathbf{L}_{1,2} \times (\mathbf{L}_{2,1} \times \mathbf{L}_{1,2})}{|\mathbf{L}_{1,2} \times (\mathbf{L}_{2,1} \times \mathbf{L}_{1,2})|}\) – единичные векторы расчетных нормалей к огибающей окружности радиуса \(R\).

Для решения полученной системы дифференциальных уравнений используется метод IDA [38].

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

Математическая модель движителей на основе гребных колес

Модель гребного колеса строится на основе наших подробных нестационарных гидро-аэродинамических расчетов методом конечных элементов с помощью программ ANSYS CFX (рис. 2), по результатам которых определялись массивы поправочных коэффициентов в аналитическую модель, аналогично методике, изложенной в [39].

Эта модель строится на основе представления тормозящего момента как

\[M_{i} = - (R\omega - v)(\rho \delta L v)R\]

\[\delta = \frac{2 (R \cos\phi + (R-\Delta H) \phi) - (R-\Delta H) \pi}{\pi - 2 \phi}\]

\[\phi = \arcsin\frac{R-\Delta H}{R}\]

где \(\phi\) – угол между поверхностью воды и радиусом в точке входа лопатки в воду, \(v\) – скорость центра колеса относительно поверхности воды в плоскости колеса, \(R\) – радиус, \(\omega\) – угловая скорость вращения, \(\Delta H\) – глубина погружения, \(L\) – ширина колеса, \(\delta\) – средняя по времени глубина погружения лопатки под невозмущенную воду.

Фактический тормозящий момент определяется из данных расчета по МКЭ исходя из соотношений \[M_{r} = \mu_{r} M_{i},\] \[\mu_{r} = \frac{M_{f}}{M_{i}^*},\] где \(M_{f}\) – момент, полученный по МКЭ, \(M_{i}^*\) – идеальный момент \(M_{i}\), рассчитанный для тех же параметров, что и \(M_{f}\).

Тогда сила в направлении движения \[F_x = \eta M_{r} \frac{\omega}{v},\] где \(\eta\) – КПД колеса, определенный по МКЭ.

Вертикальная составляющая силы колесного движителя определялась из соотношения \[F_y = F_x \frac{F_{yf}}{F_{xf}},\] где \(F_{xf,yf}\) – компоненты \(x,y\) силы, полученные по МКЭ.

\(M_f\), \(\eta\), \(F_{xf,yf}\) – функции скорости вращения колеса \(\omega\), глубины погружения \(\Delta H\), и скорости движения судна \(v\). Значения этих функций определяются интерполяцией результатов, полученных по расчету МКЭ.

Эффективность движителя \(\eta\) определялась как усредненное по периоду между входом лопаток в воду отношение мощности, связанной с силой от колеса \(\mathbf{F}\), к механической мощности на валу,

\[\eta = \frac{\mathbf{F} \cdot \mathbf{v}_{v}}{|\mathbf{M} \cdot \mathbf{\omega}|},\]

где \(\mathbf{M}\) – вращательный момент на валу, \(\mathbf{\omega}\) − угловая скорость вращения колеса.

Сила \(\mathbf{F}\) вычислялась как интеграл горизонтальных проекций сил давления и поверхностного трения по всей площади поверхности колеса и обтекателя.

Предварительно проводилась расчетная оптимизация формы лопаток и кожуха, позволившая поднять η при скорости судна \(10\ldots45\;\text{км/ч}\) до \(45\ldots 50\;\text{%}\).

Двигатель и динамика

Расчеты двигателя проводятся исходя из энергетических соображений и литературных сведений о КПД [40]. Мощность двигателя пропорциональна КПД, расходу углеводородного горючего и его теплотворной способности по [41]. Для двигателей известны зависимости КПД от частоты вращение ротора. Она рассчитывается из уравнений динамики для вращения твердого тела вокруг неподвижной оси [42].

Динамика твердого тела определяется системой уравнений [43].

\[\mathbf{v} = \mathbf{P}/M\] \[\mathbf{\omega} = \hat{I}^{-1}_{lrf}\mathbf{L}\] \[\frac{\partial \mathbf{x}}{\partial t} = \mathbf{v}\] \[\frac{\partial \mathbf{q}}{\partial t} = \frac{1}{2}\omega \mathbf{q}\] \[\frac{\partial \mathbf{P}}{\partial t} = \mathbf{F}\] \[\frac{\partial \mathbf{L}}{\partial t} = \mathbf{T}\]

Здесь \(\mathbf{P}\) – импульс тела, \(\mathbf{L}\) – момент импульса тела, \(\mathbf{F}\) – сила, действующая на тело, \(\mathbf{T}\) – момент силы, действующий на тело, \(\mathbf{q}\) – кватернион вращения тела, \(\mathbf{\omega}\) – угловая скорость вращения тела (в виде вектора), \(\mathbf{v}\) – скорость движения тела, \(M\) – масса тела, \(\mathbf{x}\) – вектор положения центра масс тела, \(\hat{I}^{-1}_{lrf}\) – матрица, обратная матрице момента инерции в СК, связанной с землей, вычисляется по формуле \[\hat{I}^{-1}_{lrf}=\hat{R}\hat{I}^{-1}\hat{R}^{T},\] где \(\hat{R}\) – матрица вращения, по известным формулам построенная из кватерниона \(\mathbf{q}\) , \(\hat{I}\) – матрица момента инерции тела в системе координат главных осей инерции тела.

При помощи вспомогательного программного обеспечения [44], полученные системы уравнений для каждой из отдельных моделей из раздела переписываются в общих терминах и упорядочиваются для вычислений. Шаг по времени для дифференциальных уравнений рассчитывается с Ньютоновским решателем для плотных матриц и шагом по времени по формуле обратного дифференцирования [45].

Математическая модель СВП с таким количеством взаимосвязанных рабочих и сопутствующих процессов разработана, по-видимому, впервые.

Надежность описанной общей комплексной модели СВП проверялась сопоставлением характеристик реального СВП классического типа [10] с соответствующей комплексной моделью, а также расчетов и экспериментов с малоразмерной физической моделью.

В интересах разработки СВП были проведены прогнозные вычисления по комплексным математическим моделям для СВП разных масштабов с классическими и инновационными ГО и движителями, различными компоновками и значениями проектируемых параметров, с разными траекториями и ускорениями над твердыми и жидкими поверхностями с разными неровностями, с препятствиями разных форм и размеров.

Экспериментальные и конструкторские разработки СВП новых типов

С помощью вышеописанных математических моделей проведен ряд экспериментальных и конструкторских разработок по программе создания СВП (рис. 3-5). Разработанные технические решения реализованы в полномасштабном опытном образце ТСВП массой около \(4,9\;\text{т}\) с гибридным опорным ходовым комплексом.

Для примера приведем некоторые технические характеристики транспортного СВП средней грузоподъемности, разрабатываемого нами в рамках ОКР «Север» (рис. 5): полная масса \(32,0\;\text{т}\), грузоподъемность \(10,0\;\text{т}\), длина \(25,9\;\text{м}\), ширина \(13,6\;\text{м}\), погрузочная высота \(1,2\;\text{м}\), площадь грузовой платформы \(127,0\;\text{м}^2\), высота преодолеваемого уступа (волны) \(2,0\;\text{м}\), преодолеваемый подъем \(21,0^\text{o}\), максимальная скорость \(45\;\text{км/ч}\).

Максимальное давление ТСВП на грунт со стороны ВП в режиме частичной разгрузки \(206\;\text{даН/м}^2\) (при поднятых колесах \(340\;\text{даН/м}^2\))

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

Иллюстрации

Литература