ВведениеПроектирование современных инженерных систем в судостроении невозможно представить без систем автоматизированного проектирования. Требования к повышению производительности труда и качества выпускаемой продукции обусловливают необходимость максимальной автоматизации всех проектных процедур, начиная с этапа технического предложения и эскиза вплоть до составления рабочей конструкторской документации.Определенную сложность для современных средств автоматизации инженерных расчетов представляют вычисления, где математическая постановка проектной процедуры неочевидна; соответственно, и последующая их алгоритмическая реализация часто неудовлетворительна [1]. Одним из таких расчетов является определение собственной частоты поперечных колебаний системы судового валопровода. Несмотря на большое количество научных работ, посвященных этой теме, даже вычисление некоторых параметров такой динамической системы представляет собой серьезную задачу. Основными факторами, влияющими на характер и интенсивность колебаний судовых валопроводов, являются статические, динамические и циклические внешние нагрузки [2]. К таким факторам относятся переменные гидродинамические изгибающие моменты, возникающие в результате работы гребного винта в неравномерном поле скоростей [3, 4]. Подобные моменты действуют в вертикальной и горизонтальной плоскостях, воспринимаются валопроводом и дейдвудным устройством и оказывают непосредственное влияние на возникновение резонансного состояния системы.В настоящей работе предложена программная реализация полуэмпирической методики определения гидродинамических изгибающих моментов, действующих на валопроводы судов. Расчет изгибающих моментов является одной из основ повышения достоверности расчетов поперечных колебаний, а использование программного пакета компьютерной алгебры Maple позволяет повысить точность и существенно сократить затраты времени инженера-проектировщика.  Материалы и методы исследованияВ работе [5] предложена полуэмпирическая методика определения гидродинамических моментов, возникающих на гребном винте при работе в неравномерном поле скоростей за корпусом судна.Определение постоянной составляющей изгибающего момента основано на результатах тензометрирования валопроводов крупнотоннажных транспортных судов в натурных условиях. Периодическая составляющая установлена на основе анализа результатов измерений полей скоростей за корпусом транспортных судов на моделях в опытовом бассейне.Исходные данные по корпусу судна для расчета, а также элементы гребного винта и величины, характеризующие режим его работы, представлены в табл. 1.Таблица 1Параметры для определения гидродинамических изгибающих моментов ПараметрОбозначениеЕдиница измеренияКоэффициент общей полнотыδ–Коэффициент формы кормовых шпангоутовτ–Диаметр гребного винтаDмЧисло лопастейz–Шаг винтаHмШаговое отношениеH/D–Дисковое отношениеA/Ad–Скорость суднаVsузVм/сЧисло оборотов винтаnmоб/минnоб/сКоэффициент попутного потокаW–Коэффициент упора винтаK1–Коэффициент крутящего момента винтаK2–Расстояние между сечением гребного вала, в котором определяется изгибающий момент от плоскости диска винтаLмСила упора гребного винта на расчетном режимеPНКрутящий момент (момент силы) гребного винта на расчетном режимеMкрН∙мПриведенная характеристика осевой неравномерностиполя скоростейα–Приведенная характеристика окружной неравномерностиполя скоростейβ–Косинусные составляющие гармоник порядка Z ± 1 в разложении поля скоростей в ряд ФурьеaZ±1–Значения тангенсов углов наклона кривых действия гребного винта в свободной воде на расчетном режиме работы винтаC1, C2–Относительная поступь гребного винта по скорости судна V/(nD)λ– Коэффициент формы кормовых шпангоутов τ может быть установлен согласно принятым в транспортном судостроении формам кормовых шпангоутов: для V-образной τ ≤ 0,1; умеренно V-образной τ = 0,1÷0,3; умеренно U-образной τ = 0,3÷0,4; U-образной τ ≥ 0,4. Также коэффициент τ может быть определен как отношение a/b (рис. 1).   Рис. 1. Построения для определения коэффициента формы кормовых шпангоутов Заметим, что вспомогательные величины C1 и C2 получены как функции шагового и дискового отношения гребного винта в результате обработки диаграмм натурных испытаний гребных винтов [6] и представлены графически на рис. 2, 3.  C1H/DH/DC2   Рис. 2. Графики для определения коэффициентов C1 и C2   Рис. 3. Графики для определения коэффициентов α и β Постоянная составляющая изгибающего момента в вертикальной плоскости:                                                                  (1)Постоянная составляющая изгибающего момента в горизонтальной плоскости:                                                       (2)Амплитуда колебаний изгибающего момента на лопастной частоте в вертикальной плоскости:                                        (3)Амплитуда колебаний изгибающего момента на лопастной частоте в горизонтальной плоскости:                                        (4)где для определения коэффициента az±1 использованы результаты анализа имеющихся экспериментальных материалов по неравномерности поля скоростей судов [7, 8], представленные в табл. 2.Таблица 2Определение коэффициента  n = z ± 1234567an–0,1800,025–0,0850,015–0,0350,005 Изгибающие моменты, действующие в вертикальной (Мво) и в горизонтальной (Мго) плоскостях в сечении гребного вала могут быть найдены из соотношений (1) и (2), а по формулам (3) и (4) с использованием приведенных в табл. 2 данных возможно выполнение расчета амплитудных значений колебаний изгибающего момента в вертикальной (Мвz) и горизонтальной (Мгz) плоскости в гребных валах для судов.  Программная реализация методики в среде Maple Произведем описание программной реализации полуэмпирической методики определения гидродинамических изгибающих моментов. Расчет изгибающих моментов выполним с помощью программного пакета компьютерной алгебры Maple. Для вычисления коэффициентов C1, C2, τ применен метод билинейной интерполяции с целью обобщения интерполяции одной переменной для функции двух переменных. Суть метода состоит в следующем. Пусть известны значения функции f в вершинах прямоугольника A11, A12, A21, A22, где точка Aij имеет координаты (xi, yj), а i, j = 1, 2 (рис. 4).   Рис. 4. Интерполированное значение исходной функции в точке (x, y) Требуется найти приближенное значение функции f в точке, лежащей внутри прямоугольника, с координатами (x, y). Для этого используем формулу      где fij = (xi, yj) – значения функции f в точке (xi, yj).Последовательность расчета приведена в табл. 3.Таблица 3Последовательность вычисления гидродинамических изгибающих моментов в среде Maple№ Этапы Программная реализация1Определениефункций согласно уравнениям (1)–(4)with(Linear Algebra):U1:=P*D*(lambda/(2*K1))*(0.355*C1*alpha+1.492*((1-W)*C2+(z*K2/lambda))*L*beta);U2:=P*D*(lambda/(2*K1))*(1.492*C2*L*alpha-0.335*((1-W)*C1+(z*K1/lambda))*beta);U3:= P*D*(lambda/(2*K1))*sqrt((0.335*C1*(a2+a1))^2+(1.492*C2*L*(a2-a1))^2);U4:= P*D*(lambda/(2*K1))*sqrt((0.335*C1*(a2-a1))^2+(1.492*C2*L*(a2+a1))^2);2Задание значений параметровK1, K2, λ, τ, L, z, W, P, D, aZ+1, aZ–1K1:=0.198:     K2:=0.029:     lambda:=0.795:     tau:=0:     L:=0.3:     z:=4:W:=0.34:     a1:=0.015:     a2:=0.025:     P:=478:    D:=4.5:3Вычисление коэффициентовC1, C2, τ методом билинейной интерполяции 3.1Задание вершинпрямоугольника матрицей АA:=Matrix(2,2,[[0.6,0.85],[0.8,1.0]]);3.2Задание значений функциив вершинах прямоугольникаматрицей FF:=Matrix(2,2,[[0.46,0.495],[0.433,0.465]]);3.3Задание интерполирующей функцииf:=(x,y)->(F[1,1]*(A[2,1]-x)*(A[2,2]-y)+F[2,1]*(x-A[1,1])*(A[2,2]-y)+F[1,2]*(A[2,1]-x)*(y-A[1,2])+F[2,2]*(x-A[1,1])*(y-A[1,2]))/((A[2,1]-A[1,1])*(A[2,2]-A[1,2]));3.4Вычисление приближенного значения функции в точкес заданными координатами С1:=f(0.7,0.9);функции C2(H/D) и α(τ) интерполируются аналогично4Задание эмпирическойзависимости β(τ) в видеинтерполяционного полиномаX:=ExcelTools:-Import("/Users/Kushner/Desktop/beta.xls", "Sheet1", "A4:A10");Y:=ExcelTools:-Import("/Users/Kushner/Desktop/beta.xls", "Sheet1", "B4:B10");for i from 1 by 1 to 7 do A:=[op(A),X[i,1]] end do:for i from 1 by 1 to 7 do B:=[op(B),Y[i,1]] end do:f:=x->interp(A,B,x);beta:=f(tau);5Задание точности вычисленийDigits:= 5:6Подстановка численныхпараметров в уравнения(1)–(4) и вычислениезначений изгибающихмоментовeval(U1);     eval(U2);     eval(U3);     eval(U4); Результат выполненного численного эксперимента по вычислению постоянных и амплитудных значений гидродинамических изгибающих моментов для судов при L = 0,3 представлен в табл. 4.Таблица 4Результаты вычислений гидродинамических изгибающих моментов Судно проектаPD, кН∙мMво, кН∙мMго, кН∙мMвz, кН∙мMгz, кН∙м1351843–25,29033,72017,7165,058134901 508–76,90854,28848,25616,588503322–5,60213,2026,7921,28815902 278–97,95445,56024,3746,83415972 154–79,69866,34323,9096,462 Рассмотренная в настоящей работе полуэмпирическая методика определения гидродинамических изгибающих моментов, возникающих на гребном винте при работе в неравномерном поле скоростей за корпусом судна, предназначена для одновальных морских судов с наиболее употребительной традиционной формой кормовой оконечности в диапазоне измерения коэффициента общей полноты корпуса δ = 0,6÷0,8 и коэффициента формы кормовых шпангоутов τ = 0÷0,6. ЗаключениеПредложенная программная реализация методики позволила алгоритмизировать и автоматизировать процесс определения внешних сил в расчетах циклической прочности гребного вала на ранних стадиях проектирования судна в условиях отсутствия данных по измерению полей скоростей в месте расположения гребного винта. Разработанная последовательность вычисления может быть применима для валопроводов многовальных судовых установок при наличии эмпирически полученных величин, характеризующих работу гребного винта на расчетном режиме. Результаты исследования являются дополнением к существующим средствам автоматизации инженерных расчетов и позволяют более подробно моделировать напряженно-деформированное состояние системы валопровода судна.