Математическое моделирование эффективных и экологичных амфибийных судов на воздушной подушке для Арктики и Дальнего Востока

Чувашев С.Н., Якимов Н.М., Чувашева Е.С., Попов С.Д., Попов С.А.

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

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

Природно-климатические и экономические особенности Арктических регионов РФ определяют особый подход к формированию транспортной системы Арктики.

Среди основных влияющих факторов:

  • низкая плотность населения (до 0,3 чел./км2), большие расстояния между малыми населёнными пунктами, что определяет относительно небольшие пассажиропотоки.

  • преимущественно очаговая деятельность при освоении сырьевых месторождений;

  • суровые климатические условия (годовой ход температур наружного воздуха в отдельных регионах достигает 100°), короткое лето, полярная ночь;

  • наличие вечной мерзлоты, значительного количества болот, рек, озер;

  • экологическая уязвимость северных экосистем;

  • повышенный уровень цен на все виды топлива, оборудование и товары в результате разового завоза из центра страны и хранения межсезонных запасов;

  • более высокая (в 7…10 раз) стоимость создания и эксплуатации традиционной транспортной инфраструктуры;

  • слабое развитие наземных транспортных коммуникаций.

Освоение минерально-сырьевых и природных ресурсов Арктики, Сибири и Дальнего Востока настоятельно требует создания адекватной транспортной системы.

Для большинства осваиваемых месторождений характерен типовой жизненный цикл «разведка – освоение – эксплуатация – свертывание» длительностью от нескольких лет до нескольких десятков лет. Как правило, добывающие компании не заинтересованы в развитии капитальных транспортных сетей, и предпочитают использовать экологически неблагополучные, малоэффективные и дорогие в эксплуатации гусеничные средства ТС, а также вертолеты, предъявляющие минимальные требования к транспортной инфраструктуре. В результате резко увеличивается стоимость добычи сырья, легко уязвимые тундровые и таежные экологические системы получают непосильную перегрузку.

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

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

Особенностью Севера является наличие многолетнемерзлых грунтов практически на всей территории. Эти грунты содержат 60…90% подземных льдов в верхнем 20…30-метровом слое. Протаивание льдов под действием естественных и антропогенных факторов приводит к образованию термокарстовых озер, бугров, провалов и передвижению почв.

Одним из самых распространенных антропогенных воздействий на природу российского Севера является нарушение ТС растительного покрова в летний период. По некоторым данным, под воздействием традиционных внедорожных ТС происходит не менее 70% всех антропогенных нарушений почвенно-растительного покрова тундровых биогеоценозов. Несущая способность переувлажненной тундры в [1] оценивалась в 10–15 кПа, что значительно меньше давления гусеничных (40…50 кПа) и особенно колёсных транспортных средств.(200…300 кПа). За несколько суток колея превращается в ручей, который спустя 2…3 года становится термокарстово-эрозионным оврагом. Естественное восстановление разрушенного почвенно-растительного покрова происходит в течение десятков, а возможно и сотен лет. В настоящее время площадь поврежденного покрова российской тундры составляет около 16%. По имеющимся оценкам, затраты на рекультивацию 1 га поврежденного оленьего пастбища составляют 30…40 тыс. руб. в ценах 1990 г. Только в Ямало-Ненецком автономном округе в результате применения транспортной техники, не адаптированной к условиям тундры, утрачено более 1 млн. 200 га пастбищ домашних оленей, а экономический ущерб природе оценивается в 60 млрд. руб. в ценах 1990 г.

Быстрый и эффективный прорыв в совершенствовании транспортного обслуживания Арктики, Сибири и Дальнего Востока может быть достигнут за счет дополнения существующей транспортной системы амфибийными средствами на воздушной подушке АСВП [2].

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

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

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

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

  1. грузоподъемность базовой машины – не менее 10 т. (предпочтительны приспособленность к размещению стандартного контейнера ISO, а также возможность сквозной погрузки транспортных средств на грузовую платформу с помощью откидных аппарелей);

  2. способность двигаться по всем видам подстилающих поверхностей со скоростью не менее 40…45 км/ч;

  3. габаритные размеры, обеспечивающие беспрепятственное движение по рекам низших технических категорий;

  4. способность уверенно и устойчиво передвигаться над неровными подстилающими поверхностями (с высотой или глубиной неровностей до 1,6…1,8 м, причем коротких или близких по длине к длине аппарата) или над водными поверхностями ( высотой волны до 1,8 м);

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

  6. способность преодолевать затяжные подъемы крутизной до 15 град., а также устойчиво двигаться по склонам крутизной до 10 град.;

  7. способность к сохранению управляемости и курсовой устойчивости при ветре со скоростью до 10…15 м/с под всеми курсовыми углами.

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

Движение традиционных АСВП обеспечивается, как правило, с помощью аэродинамического винта. Он позволяет развивать сравнительно высокие скорости (до 100 км/ч) над ровными и пологими опорными поверхностями, в т.ч. над водой, ледяной или снежной равниной. Однако, энергетическая эффективность такого движителя (отношение работы силы тяги и мощности двигателя), достаточно высокая для самолётов, для относительно тихоходных АСВП невелика, особенно резко снижается она при движении с малыми скоростями. АСВП с воздушным движителем неустойчиво и не имеет траекторной управляемости, необходимой для движения по суше, извилистым узким рекам и в иных случаях стесненной обстановки. Такие АСВП имеют крайне ограниченные возможности движения по заданной траектории даже на сравнительно небольших пологих уклонах (более 5 град).В полном объеме эти трудности можно преодолеть в амфибийном транспортном средстве с гибридным опорно-ходовым комплексом (с частичной разгрузкой контактного движителя при помощи воздушной подушки) [3–10].

При движении по воде или заболоченной местности колёсные движители малопригодны; предпочтительными представляются гребные колёса. Наши исследования показали принципиальную возможность получения высоких значений энергетической эффективности гребных колёс не только при малых скоростях (порядка 10 км/ч), но и до 40…45 км/ч [11; 12].

Воздушная подушка ВП (повышенное давление воздуха, поддерживающее аппарат) создаётся под днищем АСВП в объёме, ограниченном гибким ограждением ГО, с помощью вентиляторов-нагнетателей, приводимых в движение двигательной установкой. Гибкое ограждение классического типа расположено по периметру днища. Часто имеются внутренние перегородки, разделяющие объём воздушной подушки на 3…4 части, что повышает устойчивость АСВП. Такая конструкция ГО, правда, сильно затрудняет ремонт указанных перегородок в полевых условиях, т.к. при неработающих нагнетателях АСВП днищем опирается о землю.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Здесь c – скорость звука в воздухе, P – давление в секции подушки, V – объем секции, Gc – массовый расход на выходе нагнетателя, Qi – объемный расход в i-том элементе (всего N элементов), ρa – плотность воздуха. В рассматриваемом варианте подушка снабжена четырьмя нагнетателями, по одному на каждую секцию.

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

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

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

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

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

Нерастяжимую мнущуюся поверхность представим в виде регулярной сетки, близкой к прямоугольной. В каждом ее узле сосредоточена определенная масса (в общем случае, массы узлов могут различаться). Узлы соединены между собой невесомыми пружинами (ребрами). Искомая коническая поверхность имеет Nz+1 слоев, каждый из которых состоит из Nφ узлов. Положение и скорость (i, j)-ого узла определяется векторами xi,j и vi,j.

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

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

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

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

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

Сила трения о подстилающую поверхность \[\mathbf{F}_{fr} = -{k}_{fr}F_{r}\mathbf{\tau},\]

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

Сила давления

\[\mathbf{F} = P S \mathbf{n},\] где n – единичный вектор внешней нормали к поверхности оболочки, P – избыточное давление под оболочкой, S – площадь ячейки.

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

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

где E - модуль Юнга, например, для резины Е ≈ 107 Па; J - момент инерции сечения в соответствующем направлении. Для плоского сечения \(J = L h^3/12\), где h - толщина гибкой оболочки, L - средняя длина двух противолежащих ребер (в горизонтальном или вертикальном направлении), исходящих из данного узла, C – задаваемая координатами узлов шаблона кривизна поверхности в данной точке (в горизонтальном или вертикальном направлении), C0 – равновесные значения кривизны, при которых силы не возникают.

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

Силы для узлов в конце L1,2 рассчитываются как

\[ \mathbf{F}_{1,2} = -\frac{M}{L_{1,2}} \mathbf{n}_{1,2}, \]

а в среднем узле сила вычисляется исходя из третьего закона Ньютона:

\[\mathbf{F} = -(\mathbf{F}_1+\mathbf{F}_2).\]

Здесь n1,2 – единичные векторы расчетных нормалей к огибающей окружности радиуса R.

\[n_{1,2} = \frac{L_{1,2} \times (L_{2,1} \times L_{1,2})}{|L_{1,2} \times (L_{2,1} \times L_{1,2})|}.\]

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

Описанная методика и программа тестировались путём проверки выполнения локальных и интегральных балансовых соотношений для конического элемента ГО и его частей.

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

Рисунок 1: Наезд гибкого ограждения на полусферическое препятствие высотой 1 м со скоростью 5 км/ч (а-в) и 45 км/ч (г-е). а – наезд, вид справа, б – наезд, вид сверху, в – сход с препятствия, вид справа, г – наезд, вид сбоку-спереди-справа, д – наезд, вид сверху – спереди-справа, е – сход с препятствия, вид справа

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

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

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

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

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

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

Здесь ρa – плотность воздуха, Qi – объемный расход воздуха из i-го элемента, τi – ширина зазора под элементом, Li – длина нижней кромки данного элемента, ni – нормаль к огибающей ГО в середине элемента.

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

\[ \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} \]

Здесь απ – коэффициент заполнения площади подушки, GL – продольный коэффициент загрузки, Pavg – среднее давление в подушке, ρw – плотность воды, g – ускорение свободного падения, L – длина подушки, B – ширина подушки, S – площадь подушки, M – масса судна, κw1 – эмпирический коэффициент, задаваемый графически [20].

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

\[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}.\]

где ρw – плотность воды, v – скорость судна, τ – средняя величина зазора, L – длина подушки, S – площадь подушки, cГО=15.4 10-6, GL – продольный коэффициент загрузки, g – ускорение свободного падения, Rwave – волновое сопротивление.

Нагнетатели

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

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

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

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

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

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

\[\begin{aligned} \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}} \\ \end{aligned}\]

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

По найденным характеристикам потока, набегающего на рабочее колесо, определяются модуль и направление относительной скорости w1, число Маха набегающего потока во вращающихся координатах рабочего колеса Mw1, угол входа потока в каналы рабочего колеса β11, температура торможения Tw1*:

\[\begin{aligned} 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}{2 C_p} \\ \end{aligned}\]

где ω – угловая скорость, Cp – удельная теплоёмкость газа.

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

\[\begin{aligned} F_{21} &=& F_0 \sin \beta_{11} \\ F_{2cr} &=& \frac{{F_{21} M_{w1} ({k+1})}^{\frac{k+1}{2(k-1)}}}{\left(2+{M_{w1}^2(k-1)}\right)^{\frac{k+1}{2(k-1)}}} \end{aligned}\]

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

\[\begin{aligned} 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{2+(k-1)}{2+(k-1)M_{w2}^2}\\ a_2 &=& \frac{kk_BT_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} \\ \end{aligned}\] где kB – постоянная Больцмана, mmol – масса молекулы газа, индекс cr относится к параметрам в эффективном критическом сечении.

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

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

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

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

Значение YP находится из соотношения

\[ Y_P = \mu \chi \left(Y_{Pr} + \left(\frac{\alpha_{bin}}{\alpha_{bout}} \right)^2(Y_{Pi}-Y_{Pr})\right)\left(\frac{5t}{C}\right)^{\frac{\alpha_{bin}}{\alpha_{bout}}} + Y_S\] где значения YPr и YPi снимаются с графиков [30] в зависимости от угла лопатки на выходе αbout и от отношения шага между лопатками s к длине лопатки С, здесь αbin - угол лопатки на входе потока, t - толщина лопатки, обычно составляющая t=(0,15..0,25) С.

Влияние вязкости учитывается фактором χ, который зависит от числа Рейнольдса

\[ Re = w_2 \frac{C}{\nu} \]

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

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

при \(Re > 10^6\), \(\chi = (Re\cdot10^{-6})^{-0,2}\)

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

\[\begin{aligned} \mu &=& 1-\frac{Z_{TE}}{h} \\ Z_{TE} &=& \left(\frac{0,1 \cdot F_\theta^{0,79}}{C_R^{1/2}(h/C)^{0,55}}+32,7 \frac{\delta^*}{h}\right)h \\ \end{aligned}\]

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

\[\begin{aligned} F_0 &=& \frac{2 s}{C_x}\cos^2(\alpha_m)(\tan(\alpha_{in})+\tan(\alpha_{bout})) \\ \alpha_m &=& \frac{\alpha_{bin}+\alpha_{bout}}{2}\\ \end{aligned}\]

здесь αin - угол входа потока, Cx.- длина проекции лопатки на ось.

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

при x” от 0 до 0,3

\[ Y_S = Y_{S0} (\exp(0,9 x'')+13 x''^2+400 x''^4) \]

при x” от -0,4 до 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} \]

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

при \(x \geq 0\)

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

при \(x \leq 0\)

\[ \Delta \phi_P^2 = 1,358\cdot10^{-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 \alpha_{bin}\) – расчетный угол входа потока, din - диаметр закругления лопатки на входе.

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

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

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

Тангенциальная u1(Γ) и осевая v1(Γ) составляющие индуктивной скорости в плоскости винта представляют собой зависимости от циркуляции:

\[\begin{aligned} 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}\\ \end{aligned}\]

Затем находится скорость набегающего потока

\[W_1(\Gamma) = \sqrt{(\omega r - u_1(\Gamma))^2+(V+v_1(\Gamma))^2},\] и рассчитывается циркуляция \[\Gamma = \frac{1}{2} b C_y W_1(\Gamma).\]

Полученное уравнение решается методом бисекции относительно Γ. Расчёт начинается с концевых сечений винта. Здесь k – число лопастей, r – текущий радиус, Γ – циркуляция на данном радиусе, V – скорость осевой обдувки винта, ω – угловая скорость вращения винта, R – наибольшее значение радиуса винта, b – длина хорды лопасти. Интеграл под корнем в выражении для v1 соответствует центробежным силам. Ими, при наличии кольцевой насадки вокруг винта, можно пренебречь, т.к. она ограничивает движение в радиальном направлении. По этой же причине не учитываются концевые потери.

За винтом осевая и тангенциальная составляющие скорости:

\[\begin{aligned} u_2 &=& 2 u_1\\ v_2 &=& 2 v_1\\ \end{aligned}\]

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

Аэродинамические коэффициенты Cy, Cx определяются профилем винта. Они известны для ряда профилей из подробных экспериментальных данных [43; 44]. Эти величины находятся с помощью линейной интерполяции и экстраполяции по табличным данным для профиля. Если полагать плотность воздуха постоянной в окрестностях винта и равной ρ, подъемная сила, действующая на элемент лопасти и направленная перпендикулярно скорости набегающего потока, может быть найдена, как \(dR_1 = \rho\Gamma W_1 dr\), а сила профильного сопротивления, направленная вдоль скорости набегающего потока, как \(dR_p=\frac{C_x}{C_y} dR_1\).

Сила тяги dT1 и сила сопротивления dQ1 для элемента лопасти могут быть найдены, как проекции суммарной силы dR1+dRp на ось вращения и плоскость вращения соответственно.

Мощность dL, необходимая для вращения элемента лопасти, лежащего на окружности радиуса r, может быть определена из момента силы сопротивления dMQ1= -dQ1r:

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

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

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

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

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

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

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

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

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

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

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

\[dR_1 = \frac{1}{2} C_v ρ v^2 b_v dy,\] где Cv – аэродинамический коэффициент подъёмной силы для профиля руля, bv – длина хорды руля, v – локальная скорость набегающего на руль потока.

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

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

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

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

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

Моменты сил определяются по векторным формулам \[\begin{aligned} d\mathbf{M}_{R1} &=& \mathbf{r} \times d\mathbf{R}_1,\\ d\mathbf{M}_{Rp} &=& \mathbf{r} \times d\mathbf{R}_p. \end{aligned}\]

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

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

Движители на основе гребных колёс

Гребные колёса в течение долгого времени применялись редко из-за их малой эффективности при повышенных скоростях; однако в наших работах [11; 12] на основе вычислительных экспериментов показано, как получить высокий тяговый КПД η также при высоких скоростях (η = 45…50% при скорости судна 10…45 км/ч).

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

Тормозящий момент на валу колеса

\[\begin{aligned} M_{i} &=& -(R \omega - v)(\rho \delta L v) R, \\ \phi &=& \arcsin\left(\frac{R-\Delta H}{R}\right), \\ \delta &=& \frac{2 (R \cos \phi + (R-\Delta H) \phi) - (R-\Delta H) \pi }{\pi - 2 \phi}, \end{aligned}\]

где φ – угол между поверхностью воды и радиусом в точке входа лопатки в воду, v – скорость центра колеса относительно поверхности воды в плоскости колеса, R – радиус, ω – угловая скорость вращения, ΔH – глубина погружения, L – ширина колеса, δ – средняя по времени глубина погружения лопатки под невозмущенную воду.

Фактический тормозящий момент определяется из данных расчета по МКЭ исходя из соотношений

\[\begin{aligned} M_{r} &=& μ_{r} M_{i}, \\ μ_{r} &=& M_{f}/M_{if}, \end{aligned}\] где Mf – момент, полученный по МКЭ, Mif – идеальный момент, рассчитанный для тех же параметров, что и Mf.

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

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

Mf, η, Fxf,yf – функции скорости вращения колеса ω, глубины погружения ΔH, и скорости движения судна v. Значения этих функций определяются интерполяцией результатов, полученных по расчету МКЭ.

Эффективность движителя η определялась как усреднённое по периоду между входом лопаток в воду отношение мощности, связанной с силой от колеса F, к механической мощности на валу,

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

где M – вращательный момент на валу, ω - угловая скорость вращения колеса.

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

Рисунок 2: Гребное колесо для высоких скоростей. Цветом показана объёмная доля воды

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

Расчеты двигателей проводились исходя из энергетических соображений и литературных сведений о КПД [45]. Мощность двигателя полагается пропорциональной КПД, расходу углеводородного горючего и его теплотворной способности [46]. Характерные значения КПД задаются графически. Вращение ротора рассчитывается исходя из уравнений динамики для вращения твердого тела вокруг неподвижной оси [47]:

\[\begin{aligned} P &=& \varepsilon m' \eta(\omega) \\ M &=& \frac{P}{\omega} \\ \frac{d \omega}{d t} &=& \frac{M}{J} \end{aligned}\]

Здесь P – мгновенная мощность двигателя, ε – теплотворная способность горючего, m’ – массовый расход горючего, η(ω) – мгновенный КПД двигателя, ω – угловая скорость вращения ротора, M – крутящий момент двигателя, J – приведенный момент инерции системы относительно оси вращения.

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

\[ \begin{aligned} v_i &=& P_i/M \\ \omega_i &=& \sum_j I^{-1}_{ij} L_i \\ \frac {d x_i}{d t} &=& v_i \\ \frac {d q_0}{d t} &=& - \frac 1 2 \sum_i \omega_i q_i \\ \frac {d q_1}{d t} &=& \frac 1 2 (\omega_1 q_0 - \omega_2 q_3 + \omega_3 q_2) \\ \frac {d q_2}{d t} &=& \frac 1 2 (\omega_1 q_3 + \omega_2 q_0 - \omega_3 q_1) \\ \frac {d q_3}{d t} &=& \frac 1 2 (- \omega_1 q_2 + \omega_2 q_1 + \omega_3 q_0) \\ \frac {d P_i}{d t} &=& F_i \\ \frac {d L_i}{d t} &=& T_i \\ I^{-1} &=& R I^{-1}_{orf} R^T \\ R &=& \begin{pmatrix} q_0^2+q_1^2-q_2^2-q_3^2&2q_1q_2-2q_0q_3&2q_1q_3+2q_0q_2\\ 2q_1q_2+2q_0q_3&q_0^2-q_1^2+q_2^2-q_3^2&2q_2q_3-2q_0q_1\\ 2q_1q_3-2q_0q_2&2q_2q_3+2q_0q_1&q_0^2-q_1^2-q_2^2+q_3^2\\ \end{pmatrix} \end{aligned} \] где R - матрица вращения, Iorf – момент инерции корпуса в системе координат, связанной с ним, I – момент инерции корпуса в ЛСК, v – скорость, x – координаты центра масс в ЛСК, P – импульс, M – масса судна, L – момент импульса, q – кватернион вращения, ω – угловая скорость вращения, F – сумма сил, действующих на корпус, T – сумма моментов, действующих на корпус, i,j = 1,2,3 – индексы компонент векторов.

Методика расчета комплекса

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

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

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

Результаты вычислительных экспериментов

Расчеты по вышеописанной комплексной математической модели АСВП применялись для исследования различных типов АСВП; ниже приведены результаты вычислений для примеров:

  • Легкого АСВП классического типа (АСВП КТ) массой 18 тонн (освоенный диапазон параметров АСВП – для верификации математической модели)

  • Тяжелого АСВП КТ массой 200 тонн (для математического обеспечения при разработки нового АСВП в рамках НИР «Амфибия»)

  • Инновационного АСВП с контактными колёсными движителями, гребными колёсами и ГО из отдельных конусов.

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

Волны на глубокой воде моделируются в соответствии с [51]. В приводимых далее примерах рассматривались волны для случая скорости ветра 2.7 м/с на высоте 10 м над уровнем моря. Неровности твердой поверхности задавались с помощью генератора случайных чисел с максимальным перепадом между вершинами и впадинами 20 см и характерным размером по горизонтали 1 м и 5м (малоразмерные и крупные неровности).

Рисунок 3: Форма поверхности. a – волнующейся воды в начальный момент времени, б – с мелкими неоднородностями, в – с крупными неоднородностями

Легкое АСВП КТ

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

Гибкое ограждение имеет прямоугольную форму и состоит из 60 элементов. Ширина ГО 8.2 м, длина 19 м, высота 1.7 м. Угол наклона навесных элементов θ=45⁰. Ширина навесных элементов – 30 см. В воздушном тракте за нагнетателями характерная поперечная площадь сечения принималась равной 1 м2, а характерная длина пути – 10 м.

АСВП КТ снабжено двумя винтовыми воздушными четырехлопастными движителями с радиусом 3 м. Длина хорды лопасти винта 0.2 м, радиус кока 0.35 м, угол крутки φ(r) = (0.15 R/r + π/10) рад. Каждый винт имеет по три лопатки рулей с длиной хорды 0.3 м. В данном примере выбран профиль лопаток СибНИА С-11.

Рассматривались варианты с аэродинамическим коэффициентом лобового сопротивления cx = 0.7 и 0.35.

Режим работы двигателя для нагнетателей подбирался исходя из стабильной штатной работы воздушной подушки. Для всех вариантов подстилающей поверхности, парение стабильно при расходе горючего на работу компрессоров около 72 кг/ч. Средняя ширина зазора составляет 13 см при парении над водой и 10 см при парении над твердой поверхностью.

При разгоне задавалось нарастание расхода горючего на работу движителей. Первые 10 секунд до 90% максимального значения, и далее плавный выход на максимум 54 кг/с в течение 20 секунд. При этом движители выходят на квазистацонарный режим с задержкой 7..10 с. Сила тяги проходит через максимум из-за особенностей взаимодействия винта с набегающим потоком. В начале разгона (3..5 с) наблюдаются колебания АСВП в вертикальном направлении и соответствующие колебательные изменения зазоров гибкого ограждения, связанные с отличием значений стартовых параметров от равновесных. С ростом скорости импульсное и профильное сопротивления возрастают, а остаточное, волновое и струй – проходят через максимум, но суммарное сопротивление возрастает монотонно (рис. 5). За 60 секунд судно разгоняется до скорости 34.5 км/ч (рис. 4, а). Квазистационарная скорость при данном расходе топлива составляет 45 км/ч.

Рисунок 4: Скорость легкого АСВП КТ. а-в при разгоне по ровной водной поверхности, а – в стандартном варианте, б – при пониженном расходе топлива, в – с пониженным аэродинамическим сопротивлнением; г, д – при маневрировании на ровной поверхности, г – водной, д – твердой
Рисунок 5: Сила тяги и сопротивления легкого АСВП КТ при разгоне по ровной водной поверхности. а – сила тяги, б-ж – сопротивления, б – суммарное, в – импульсное, г – профильное, д – струй, е – волновое, ж – остаточное.

При понижении максимального расхода топлива до 27 кг/ч, АСВП не преодолевает максимум волнового и остаточного сопротивления, что приводит к стаблилизации скорости на уровне 24 км/ч (рис. 4, б)

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

При маневрировании на различных подстилающих поверхностях, расход горючего двигателем движителя фиксируется на уровне 54 кг/ч. Рули управляются ПИД-регулятором , оперирующим углом рыскания судна с параметрами P=1.2, I=0.1, D=3.6.

Типичные траектории движения приведены на рис. 7

Рисунок 6: Сила тяги и сопротивления легкого АСВП КТ при разгоне по ровной водной поверхности для пониженного коэффициента аэродинамического сопротивления. а – сила тяги, б-ж – сопротивления, б – суммарное, в – импульсное, г – профильное, д – струй, е – волновое, ж – остаточное.
Рисунок 7: Траектории движения легкого АСВП КТ при маневрировании. Стреками обозначена ориентация судна. а – на гладкой воде, б – на волнующейся воде, в – на твердой поверхности с мелкими неоднородностями, г – на твердой поверхности с крупными неоднородностями

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

Рисунок 8: Сила тяги и сопротивления легкого АСВП КТ при маневрировании на ровной водной поверхности. а – сила тяги, б-ж – сопротивления, б – суммарное, в – импульсное, г – профильное, д – струй, е – волновое, ж – остаточное.

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

Маневрирование на твердой ровной поверхности отличается от ровной водной отсутствием волнового и остаточного компонент сопротивления. В результате, при маневрировании судно разгоняется до 66 км/ч (рис. 4, д)

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

Рисунок 9: Крен, тангаж, рыскание судна. а-в – легкого АСВП КТ при маневрировании над твердой поверхностью с мелкими неоднородностями, г-е – тяжелого АСВП КТ при разгоне над твердой поверхностью с мелкими неоднородностями, ж-и – тяжелого АСВП КТ при разгоне над волнистой водной поверхностю при поперечном ветре

Тяжелое АСВП КТ

В этом примере гибкое ограждение имеет прямоугольную форму и состоит из 200 элементов. Ширина ГО 15 м, длина 35 м, высота навесных элементов 1 м. Угол наклона навесных элементов θ=45⁰. Ширина навесных элементов – 50 см.

АСВП снабжено четырьмя винтовыми воздушными четырехлопастными движителями с радиусом 1.5 м, расположенными по углам корпуса. Длина хорды лопасти винта 0.5 м, радиус кока 0.35 м, угол крутки φ(r) = 0.08 R/r рад. Поворот осуществляется поворотом движителей вокруг вертикальной оси.

Аэродинамический коэффициент лобового сопротивления cx = 0.4, аэродинамический коэффициент поперечного сопротивления cx=0.5.

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

Средняя ширина зазора при парении над ровной поверхностью 3,5 см, над ровной водной поверхностью – 7,7 см.

Ниже приведены результаты расчетов разгона на волнующейся воде для различных направлений ветра (встречный и поперечный).

При разгоне по ровной твердой поверхности, судно разгоняется до крейсерской скорости 72 км/ч за 60 секунд (рис. 10, а). Расход топлива на работу движителей на крейсерской скорости – 72 кг/ч.

Рисунок 10: Скорость тяжелого АСВП КТ при разгоне. а,б – по ровной поверхности, а – твердой, б – водной; в – по твердой с мелкими неоднородностями (сливается с а), г – по волнистой водной при боковом ветре

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

При разгоне по ровной водной поверхности, за 60 секунд судно достигает половины крейсерской скорости и далее постепенно разгоняется до крейсерской в течение 3-4 минут (рис. 10, б). При этом вклад гидродинамических компонент сопротивления в суммарное оказывается решающим (рис. 11). Расход топлива в процессе разгона в среднем составляет 612 кг/ч. Расход топлива на крейсерской скорости – 382 кг/ч.

Рисунок 11: Cопротивления при разгоне тяжелого АСВП КТ по ровной водной поверхности. а – суммарное, б – импульсное, в – профильное, г – струй, д – волновое, е – остаточное.

Примечательно, что в данном примере суммарное сопротивление проходит через максимум при скоростях около 40 км/ч, т.е. для повышения топливной эффективности двигаться АСВП по воде следует либо с малыми скоростями – до 36 км/ч, либо быстрее 54 км/ч.

При разгоне на волнующейся воде, возникают колебания характеристик, связанные с неровностью поверхности. Это приводит к значительной тряске (рис. 9, ж, з), однако в остальном влияет на поведение судна слабо (рис. 10, г; рис. 9, и]). При поперечном ветре, возникает дрейф, связанный с вкладом профильного сопротивления ветру (поперечная скорость достигает 3,6 км/ч). При встречном ветре, возрастает вклад профильной и импульсной компонент сопротивления, и преодоление пика сопротивления оказывается затруднено. В результате, разгон до крейсерской скорости занимает 16 минут, а расход топлива на крейсерской скорости возрастает до 418 кг/ч.

Рассматриваись сценарии маневрирования над гладкими водной и твердой поверхностями, над поверхностю с крупными неоднородностями и над волнующейся водой при попутном ветре 2,5 м/с. Траектория движения для всех подстилающих поверхностей приблизительно одинакова. Типичный образец приведен на рис. 12.

Рисунок 12: Траектория движения тяжелого АСВП КТ при маневрировании на ровной твердой поверхности. Стреками обозначена ориентация судна

Видно, насколько АСВП инерционно: поворот происходит на расстоянии порядка 5 км, при поворотах скорость судна снижается мало, для выполнения поворота АСВП долгое время должно двигаться под углом к курс�� порядка 45⁰. При маневрировании над твердой поверхностью скорость судна варьируется в диапазоне от 70 до 73 км/ч, а над водой – от 68 до 69 км/ч.

В ряде случаев при более резких поворотах моделировалась аварийная ситуация – переворот АСВП.

Инновационное АСВП для российского Севера

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

На рис. 14-17 для такого инновационного АСВП представлены некоторые результаты вычислений по комплексной модели. Масса АСВП 50200 кг, аэродинамические коэффициенты по осям x и y – 0.4 и 0.5. Основные характеристики ходовой части: Ширина и диаметр гребных колес: 4 м, Момент инерции 0.0146 кг·м². Для обычных колес диаметр 4 м. Момент инерции 7200 кг·м². Коэффициент реакции поверхности 107 Н/м. Коэффициент трения поверхности 1. Коэффициент вязкого трения вращения в колесах 1 1/с. Коэффициент амортизации 2500 Н·с/м Коэффициент жесткости пружины амортизатора 250000 Н/м. Масса колеса 2000 кг. Коэффициент передачи с двигателя колеса на само колесо 25. Диаметр нижнего торца большого конического элемента ГО 6.04 м, малого – 4.2 м.

На твердой поверхности при скорости 45… 47 км/ч поворот происходит на расстоянии порядка 20 м, как у автомобиля. Бертеновские конуса создают 120 кН подъёмной силы, колёсные контактные движители – 360 кН, при этом избыточное давление в ГО составляет порядка 1 кПа, т.е. 1% от атмосферного. САУ задвижки существенно сокращает колебания ГО и АСВП в целом, вертикальные колебания амплитудой 5 мм, частотой 1 Гц. Расход горючего двигателей колёсных контактных движителей составляет 230,4 кг/ч.

Хорошая управляемость такого АСВП продемонстрирована также на других траекториях и формах твёрдых поверхностей, в т.ч. гладких с уклонами, с малоразмерными и крупными неоднородностями.

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

Рисунок 13: Пример схемы ходовой части инновационного АСВП. 1 – конические элементы ГО, 2 – контактные колесные движители, 3 – гребные колеса
Рисунок 14: Скорость инновационного АСВП при маневрировании на твердой поверхности
Рисунок 15: Траектория движения инновационного АСВП при маневрировании на твердой поверхности
Рисунок 16: Мощность на валах колёс. а – переднего правого, б – переднего левого, в – заднего правого, г – заднего левого
Рисунок 17: Давление в переднем правом элементе ГО

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

Некоторые результаты расчётов движения АСВП по воде показаны на рис. 18 для примера разгона по гладкой воде от средних скоростей к высоким. Видно, что АСВП разгоняется за несколько секунд, средний упор гребных колёс плавно растет с увеличением скорости и связанных с этим сопротивлений, при этом КПД примерно постоянен (около 50%). Несмотря на значительные выбросы во временных зависимостях сил при гребках, центр масс вибрирует незначительно.

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

Рисунок 18: Временные зависимости вертикальной координаты центра масс АСВП (a), скорости АСВП (б), упора колёс в горизонтальном направлении (в) и мощности двигателя (г) для разгона инновацонного АСВП по водной проверхности

Заключение

Из анализа природных, экологических и социально-экономических условий применения транспортных средств на крайнем Севере следует перспективность разработки амфибийных судов на воздушной подушке АСВП. Разработан математический аппарат для концептуального проектирования как АСВП классического типа, так и АСВП, содержащих инновационные технические решения: гибкое ограждение, состоящее из отдельных конусов, колёсные контактные движителии новые, высокоэффективные гребные колёса с высоким тяговым КПД как на малых, так и на высоких скоростях (до 45 км/ч). Для условий Севера показана целесообразность разработки инновационного АСВП.

Литература

1. Чеботаев А.А. Безвредные Транспортные Средства Для Севера / А.А. Чеботаев, А.Д. Мельник // Сб. Матерталов Всесоюзной Научно-Практической Конф. 23–26 Октября 1990 Г Научно-технический прогресс и перспективы развития новых специализированных видов транспорта. – М.: ВНИИПК техоргнефтегазстроя, 1990. – Vol. 2. – P. 115-125.

2. Попов С.Д. Фундаментальные Проблемы Развития Внутреннего Транспорта Малонаселенных Регионов России и Пути Ее Решения На Основе Развития Высокомобильных Автомобильных Транспортных Комплексов (На Примере Архангельской Области) / С.Д. Попов // Сб.Статей По Итогам Международной Научно-Практической Конференции. Экономика, проектный менеджмент, образование, юриспруденция, экология, медицина, социология, философия, филология, психология, техника, математика: состояние и перспективы развития. – СПб: КульИнформПресс, 5BC–2013. – P. 110-117.

3. Попов С.Д. Проектирование и Комплексное Математическое Моделирование Судна На Воздушной Подушке Для Регионов Севера, Сибири и Арктического Континентального Шельфа. / С.Д. Попов, С.Н. Чувашев // Инженерный журнал: наука и инновации. – 2013. – № 3 (15). – P. 9.

4. Belousov B.N. Heavy-duty wheeled vehicles: Design, theory, calculations. Heavy-duty wheeled vehicles / B.N. Belousov. – Warrendale, Pennsylvania, USA: SAE International, 2014. – 553 p.

5. Попов С.Д. Некоторые Проблемы Создания Амфибийных Транспортных Систем, Предназначенных Для Решения Транспортных Задач и Освоения Труднодоступных Регионов Севера и Сибири, а Также На Арктическом Шельфе. / С.Д. Попов // Материали За 9-а Международна Научно-Практична Конференция Achievement of high scool. – София «Бял ГРАД-Ы» ООД, 2013. – Vol. 45. – P. 104.

6. Попов С.Д. Разработка Технологии Выбора Несущего Комплекса Для Транспортных Средств На Воздушной Подушке (ТСВП), Предназначенных Для Эксплуатации На Севере и в Сибири / С.Д. Попов, С.Н. Чувашев // Сб. Статей Международной Научно-Практической Конференции Инновационное развитие современной науки / ed. А.А. Сукиасян. – Уфа: РИЦ БашГУ, 2014. – Vol. 3. – P. 287-295.

7. Попов С.Д. Разработка Технологии Математического Моделирования Некоторых Опасных Ситуаций При Эксплуатации ТСВП, Предназначенных Для Эксплуатации На Севере и в Сибири. / С.Д. Попов, С.Н. Чувашев // Сб. Статей Международной Научно-Практической Конференции Инновационное развитие современной науки / ed. А.А. Сукиасян. – Уфа: РИЦ БашГУ, 2014. – Vol. 3. – P. 300-306/366.

8. Попов С.Д. Разработка Технологии Выбора Движительных Комплексов Повышенной Эффективности Для ТСВП, Предназначенных Для Эксплуатации На Севере и в Сибири. / С.Д. Попов, С.Н. Чувашев // Сб. Статей Международной Научно-Практической Конференции Инновационное развитие современной науки / ed. А.А. Сукиасян. – Уфа: РИЦ БашГУ, 2014. – Vol. 3. – P. 296-230/366.

9. Попов С.Д. Отработка Технологий Исследований Составных Частей и Моделей ТСВП, Предназначенных Для Эксплуатации На Севере и в Сибири. / С.Д. Попов, К.В. Долотов, Б.В. Овсянников // Сб. Статей Международной Научно-Практической Конференции Инновационное развитие современной науки / ed. А.А. Сукиасян. – Уфа: РИЦ БашГУ, 2014. – Vol. 3. – P. 281-286/366.

10. Попов С.Д. Амфибийные Транспортные Средства с Гибридным Опорно-Ходовым Комплексом. / С.Д. Попов, А.Е. Дубин // Русский инженер-транспортник (авиация, автомобили, спецтехника). – P. 36-37.

11. Якимов Н.М. О Возможности Повышения Эффективности Гребного Колеса При Высокой Скорости Судна : Машиностроение / Н.М. Якимов, С.Д. Попов, С.Н. Чувашев // Вестник Московского государственного технического университета им. Н.Э. Баумана. – 2016. – № 1. – P. 101-111.

12. Якимов Н.М. Выбор Геометрии Эффективного Высокоскоростного Судового Движителя На Основе Гребного Колеса с Учетом Гидроаэродинамики : Машиностроение / Н.М. Якимов, С.Д. Попов, С.Н. Чувашев // Herald of the Bauman Moscow State Technical University. Series Mechanical Engineering. – 2015. – № 6. – P. 50-58.

13. Чувашева Е.С. Комплексная Математическая Модель Для Концептуального Проектирования Высокоскоростных Летательных Аппаратов / Е.С. Чувашева, С.Н. Чувашев, И.Г. Зорина // Информационные технологии. – 2012. – № 11(195). – P. 10-14.

14. Чувашева Е.С. Выбор Рациональных Характеристик Высокоскоростных Летательных Аппаратов Разных Масштабов На Основе Комплексной Математической Модели / Е.С. Чувашева, С.Н. Чувашев // Информационные технологии. – 2013. – № 8. – P. 12-16.

15. Демешко Г.Ф. Проектирование Судов. Амфибийные Суда На Воздушной Подушке Книга 2 / Г.Ф. Демешко. – СПб: Судостроение, 1992. – 329 p.

16. Gravdahl J.T. Drive torque actuation in active surge control of centrifugal compressors / J.T. Gravdahl, O. Egeland, S.O. Vatland // Automatica. – 2002. – Vol. 38. – № 11. – P. 1881-1893.

17. Hindmarsh A.C. The PVODE and IDA algorithms / A.C. Hindmarsh. – Technical Report UCRL-ID-141558, LLNL, 2000.

18. Демешко Г.Ф. Проектирование Судов. Амфибийные Суда На Воздушной Подушке Книга 1 / Г.Ф. Демешко. – СПб: Судостроение, 1992. – 269 p.

19. Прохоров А.М. Физический Энциклопедический Словарь / А.М. Прохоров. – М.: Советская энциклопедия, 1983.

20. Справочник По Теории Корабля. Vol. 3. Управляемость водоизмещающих судов. / ed. Я.И. Войткунский. – Л.: Судостроение, 1985.

21. Mantle P.J. Air cushion craft development, first revision / P.J. Mantle. – DAVID W TAYLOR NAVAL SHIP RESEARCH AND DEVELOPMENT CENTER BETHESDA MD, 1980.

22. Шерстюк А.Н. Насосы, Вентиляторы, Компрессоры / А.Н. Шерстюк. – М.: Высшая школа, 1972. – 343 p.

23. Турбины Тепловых и Атомных Электрических Станций / eds. А.Г. Костюк, В.В. Фролов. – М.: МЭИ, 2001. – 488 p.

24. Boyce M.P. Gas turbine engineering handbook / M.P. Boyce. – Amsterdam ; Boston: Elsevier/Butterworth-Heinemann, 2012. – 956 p.

25. Moroz L. A uniform approach to conceptual design of axial turbine / compressor flow path / L. Moroz, Y. Govorushchenko, P. Pagur, [et al.] // Future of Gas Turbine Technology 3 rd International Conference 11-12 October 2006.

26. Kikstra J.F. Dynamic Modeling of a Cogenerating Nuclear Gas Turbine Plant—Part I: Modeling and Validation / J.F. Kikstra, A.H.M. Verkooijen // Journal of Engineering for Gas Turbines and Power. – 2002. – Vol. 124. – Dynamic Modeling of a Cogenerating Nuclear Gas Turbine Plant—Part I. – № 3. – P. 725-733.

27. Абрамович Г.Н. Прикладная Газовая Динамика. Vol. 1 / Г.Н. Абрамович. – М.: Физматгиз, 1991. – 600 p.

28. Абрамович Г.Н. Прикладная Газовая Динамика. Vol. 2 / Г.Н. Абрамович. – М.: Физматгиз, 1991. – 304 p.

29. Черный Г.Г. Газовая Динамика / Г.Г. Черный. – М.: Наука. Гл. ред. физ.-мат. лит., 1988. – 424 p.

30. Ainley D.G. An Examination of the Flow and Pressure Losses in Blade Rows of Axial-flow Turbines : ARC technical report / D.G. Ainley, G.C.R. Mathieson, A.R.C. Ministry of Supply. – H.M. Stationery Office, 1951.

31. Benner M.W. Influence of Leading-Edge Geometry on Profile Losses in Turbines at Off-Design Incidence: Experimental Results and an Improved Correlation / M.W. Benner, S.A. Sjolander, S.H. Moustapha // Journal of Turbomachinery. – 1997. – Vol. 119. – Influence of Leading-Edge Geometry on Profile Losses in Turbines at Off-Design Incidence. – № 2. – P. 193-200.

32. Александров В.Л. Воздушные Винты / В.Л. Александров. – М.: Оборонгиз, 1951.

33. Кравец А.С. Характеристики Воздушных Винтов / А.С. Кравец. – М.: Оборонгиз, 1941.

34. Ветчинкин В.П. Теория и Расчет Воздушного Гребного Винта / В.П. Ветчинкин, Н.Н. Поляков. – М.: Оборонгиз, 1940.

35. Шайдаков В.И. Аэродинамическое Проектирование Лопастей Воздушного Винта / В.И. Шайдаков, А.Д. Маслов. – М.: МАИ, 1995. – 34 p.

36. Шайдаков В.И. Аэродинамический Расчет Вертолета / В.И. Шайдаков. – М.: МАИ, 1988.

37. Юрьев Б.Н. Аэродинамический Расчет Вертолета / Б.Н. Юрьев. – М.: Оборонгиз, 1956.

38. Sparenberg J.A. Hydrodynamic propulsion and its optimization: Analytic theory : Fluid mechanics and its applications. Vol. 27. Hydrodynamic propulsion and its optimization / J.A. Sparenberg. – Dordrecht ; Boston: Kluwer Academic Publishers, 1995. – 368 p.

39. Larrabee E.E. The screw propeller / E.E. Larrabee // Scientific American. – 1980. – Vol. 243. – P. 134.

40. Hughes M.J. Analysis of Multi-component Ducted Propulsors in Unsteady Flow / M.J. Hughes. – Massachusetts Institute of Technology, 1993.

41. Ardito Marretta R.M. Hybrid numerical technique for evaluating wing aerodynamic loading with propeller interference / R.M. Ardito Marretta, G. Davi, G. Lombardi, [et al.] // Computers & Fluids. – 1999. – Vol. 28. – № 8. – P. 923-950.

42. Molland A.F. The maritime engineering reference book a guide to ship design, construction and operation / A.F. Molland. – Amsterdam; Boston; London: Butterworth-Heinemann, 2008.

43. Остославский И.В. Характеристики Винтовых Профилей Типа Кларк-У / И.В. Остославский, Д.В. Халезов // Технические заметки ЦАГИ. – 1937. – № 154.

44. Кашафутдинов С.Т. Атлас Аэродинамических Характеристик Крыловых Профилей / С.Т. Кашафутдинов, В.Н. Лушин. – Новосибирск: СибНИА, 1994. – 78 p.

45. Denery T. Multi-Domain Modeling of the Dynamics of a Hovercraft for Controller Development / T. Denery // AIAA Modeling and Simulation Technologies Conference and Exhibit : Guidance, navigation, and control and co-located conferences. – American Institute of Aeronautics and Astronautics, 2005.

46. Прохоров А.М. Физическая Энциклопедия. Vol. 5 / А.М. Прохоров. – М.: Большая Российская энциклопедия, 1998.

47. Кошкин Н.И. Механика. Динамика Вращательного Движения / Н.И. Кошкин // Элементарная Физика: Справочник. – М.: Наука. Гл. ред. физ.-мат. лит.,. – P. 32.

48. Witkin A. Physically Based Modeling: Principles and Practice – Constrained Dynamics / A. Witkin // COMPUTER GRAPHICS. – 1997. – Physically Based Modeling. – P. 11-21.

49. Якимов Н.М. Программное Средство Для Комплексного Математического Моделирования Сложных Технических Объектов / Н.М. Якимов, С.Н. Чувашев // ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ. – 2014. – № 11. – P. 23-30.

50. Cohen S.D. CVODE, a stiff/nonstiff ODE solver in C / S.D. Cohen, A.C. Hindmarsh // Computers in physics. – 1996. – Vol. 10. – № 2. – P. 138-143.

51. Pierson W.J. A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii / W.J. Pierson, L. Moskowitz // Journal of Geophysical Research. – 1964. – Vol. 69. – № 24. – P. 5181-5190.