CC BY
100
Еще раз об уравнении управляемости Номото Р.Г. Степахно
Судоводительский факультет МГТУ, кафедра управления судном и промыслового рыболовства
Аннотация. Предлагается новый подход к обработке результатов натурных экспериментов по стандартному маневру "Зигзаг". Модель, которая при этом идентифицируется, есть простейшая линейная модель, носящая название модели Номото (или Номото-Норрбина). Она широко используется в зарубежной практике обработки, хотя в отечественном судостроении и эксплуатации флота отношение к ней неоднозначное. В значительной мере это связано с неудобством традиционных приемов обработки таких экспериментальных данных. Именно этот вопрос обсуждается в настоящей работе, где демонстрируется иной подход как к получению экспериментальных данных, так и к их обработке в направлении ее упрощения.
Abstract. In the paper the new approach to processing of the outcomes of full-scale experiments on the standard manoeuvre "Zigzag" has been considered. The model which thus has been identified is an elementary linear model called the Nomoto model (or NoMOTo-Norrbin). It is widely used in the processing practice abroad, though in domestic shipbuilding and maintenance of fleet the relation to it is ambiguous. Largely it is connected to inconvenience of conventional techniques of these experimental data processing. Just this problem has been considered in the work, the author has demonstrated another approach both to obtaining of experimental data and to their processing trying to simplify it.
1. Введение
Постановка натурных экспериментов по маневрированию морских судов и обработка их результатов подчинены основной задаче - оценке параметров выбранной модели судна (Басин, 1977). Это сложная задача, особенно с учетом того, что сама модель обычно нелинейна. Поэтому исследователи стремятся упростить решение задачи, выбирая модели настолько простые, насколько это возможно без потери в отражении характерных черт поведения реального объекта. В 30-е годы XX в. это особенно тесно связывалось с попыткой создать стандартные методы испытаний управляемости судов и ввести сравнительные количественные критерии этого качества (Kempf, 1932).
Новую жизнь этой идее дал Номото (Nomoto et al, 1957; Nomoto and Norrbin, 1969), который предложил использовать простейшую модель управляемости судна, содержащую лишь два параметра. Эти параметры он находил по результатам специального испытания судна - маневра "Зигзаг". Для их получения Номото использовал методику, которая затем признавалась слишком сложной, громоздкой (Кацман и др., 1970). Впрочем, при существовавшем в те годы уровне развития вычислительной техники такая оценка вполне правомочна. Хотя последующие исследования других авторов (Асиновский, 1969) предлагали некоторые упрощения в обработке результатов и определении параметров модели, эти методы не получили должного распространения в отечественной судостроительной и эксплуатационной практике. Подвергалась критике (Асиновский, Гофман, 1967) и сама упрощенная модель Номото, параметры которой зависят не только от самого судна, но и от вида маневренного испытания. Тем не менее, интерес к модели Номото и испытанию "Зигзаг" в мире сохранился, особенно если учесть большое количество экспериментального материала, связанного с ними.
Именно в силу этого интереса рассмотрим другой метод параметрической идентификации модели Номото, используя тот же типовой маневр.
2. Аналитическое решения задачи Номото
Модель Номото имеет простейший вид линейного дифференциального уравнения 1-го порядка,
Trfa/dt + со= Щ (1)
имеющего два параметра T и K. Параметр T характеризует устойчивость судна на курсе, параметр К - его поворотливость в начальный момент ухода судна с постоянного курса. Само уравнение устанавливает дифференциальную зависимость угловой скорости о поворота судна вокруг вертикальной оси от кладки руля S. Для идентификации этих параметров используется маневр "Зигзаг", который состоит в следующем. На судне, идущем на прямом курсе, руль кладут, скажем, на 10° на правый борт. Судно начинает поворачивать вправо и когда изменение курса (в данном случае курс увеличивается) достигнет того же значения в 10°, руль перекладывают на 10° влево. Судно некоторое время продолжает поворот направо,
замедляет его, затем начинает поворот налево. Когда изменение курса по отношению к начальному (теперь курс уменьшается) становится равным 10°, руль перекладывают вправо и т.д. С некоторого момента этот процесс можно считать приближенно периодическим, если отвлечься от действия внешних условий плавания - ветра и течения. Естественно, что для выполнения такого маневра следует выбрать "время и место".
Маневр можно представить графически (рис. 1). На нем хорошо видно, как изменяется кладка руля в течение маневра и угол отклонения курса от своего первоначального значения. Требование при выполнении маневра равенства изменения курса и кладки руля выражается в том, что кривая их графического представления во времени имеют общие точки, например, А и В.
При существующих методиках обработки результатов требуется запись курса судна через каждые 510 сек. Разумеется, эти данные всегда полезны, особенно если они получены с самописца. Еще лучше, если такая запись ведется и для угловой скорости поворота. Но для прелагаемой методики достаточно фиксировать в эксперименте только два момента времени: время уравнивания изменения курса и кладки руля (а это существенно необходимо для реализации самого испытания) и время, когда судно возвращается на первоначальный курс. Если бы процесс был идеально периодическим, достаточно наблюдать эти моменты лишь один раз. Но лучше сделать это несколько раз, т.е. выполнить процедуру перекладки руля, скажем, 6-кратно. Тогда можно оценить степень непериодичности процесса и взять средние значения наблюдаемых моментов времени.
Чтобы свободно оперировать далее кладкой руля в модели (1), разложим ее как функцию времени в ряд Фурье. При этом ясно, что это будет разложение только по синусам, т.к. функция <5(г) нечетная. Кроме того, в ряде будут присутствовать только члены с нечетными номерами, поскольку на втором полупериоде (Т, 27) значения кладки получены зеркальным отражением в оси времени кладки первого полупериода (0, 7). Итак, если полупериод функции есть Т, то разложение кладки руля таково:
Рис. 1. Графическое представление во времени кладки руля, отклонения курсового угла и угловой скорости маневра "ЗИГЗАГ" (Т = 75 с, ^ = 10 с, г3 = 45 с, Т1 = 7.55 с)
¿¡(г) = Е^^Ь^т(пк4/Т), к = 1, 3, 5 Ьк = 2 [ [ ¿(0 ^т(пк4/Т)ёг]/Т.
(2)
Если теперь представить кладку руля как кусочно-заданную функцию на первом полупериоде в форме
Гйм, г < п ¿(0 = <¡¿0, л < г < т-и
[¿0 (Т-о/л, т-г1 < г < т,
то, взяв от этой функции интеграл (2), получим выражения для коэффициентов разложения Фурье кладки руля
Ьк = 4^ш(лк-/1/Т) / [(лк/Т)2Т-г1].
(3)
Полученный ряд Фурье неплохо сходится (со скоростью, обратной квадрату номера к); в подтверждение этого приведем в табл. 1 результаты восстановления значений кладки руля для 3, 6 и 9 членов в сумме ее представления рядом (2). В табл. 1 выбраны для расчета Т = 75 с, и = 10 с, ¿>0 = 10°. Представлена только левая половина полупериода (0, 35), т.к. правая половина полностью симметрична относительно средней точки Т/2 = 37.5 с интервала (0, Т).
Хорошо видно, как улучшается сходимость ряда к точному значению кладки, приведенному в нижнем ряду первой строки табл. 1 при увеличении числа членов ряда. Разумеется, от явлений колебаний Гиббса, особенно вблизи точек и Т-^ разрыва производной функции не удается избавиться полностью. Но подобная точность представления нас вполне устраивает, особенно применительно к решаемым далее задачам.
Таблица 1
г, с 0 5 10 15 20 25 30 35
N /8, град 0 5 10 10 10 10 10 10
3 0 5.37 8.99 10.32 10.16 9.80 9.87 10.15
6 0 5.10 9.40 10.18 9.85 10.12 9.92 10.03
9 0 5.02 9.56 9.98 10.06 9.94 10.04 9.98
Считая известным представление Фурье функции (2), будем искать решение уравнения (1) в виде более общего тригонометрического ряда
®(0 = Ък[Л,соэ(пк4/Т) + В^ш (пк4/Т)], (6)
а с его помощью выразить и изменение курсового угла <рна основании известного соотношения <а = с1фЛ:
р(0 = Т,к[Лк$Щ%к4/Т - В£яъ(%к4/Т)]/(%к/Т). (7)
Подставляя предполагаемое решение и его производную в уравнение (1) вместе с рядом, представляющим его правую часть, и уравнивая коэффициенты при синусах и косинусах одной частоты, найдем в итоге этой обычной процедуры коэффициенты решения Лк и Вк:
Лк = -К Ьк Т (пк/Т) / [1 + (пкТ^/Т)2] Вк = КЬк / [1 + (лкВД2].
2-|
(8)
Это позволяет записать аналитическое решение в форме единой функции, а не пользоваться кусочным представлением решения, как это делалось обычно для данной модели. Именно это обстоятельство позволяет решить задачу идентификации параметров 71 и К модели (1).
3. Идентификация параметров модели
По самому существу маневренных испытаний вида "Зигзаг" мы имеем два граничных условия -равенство изменения курсового угла и угла кладки руля (с учетом знаков) в зафиксированные при эксперименте моменты времени ^ и (Г-^). Правда, в силу цикличности процесса, оба эти граничных условия приводят к одному и тому же уравнению, в которое входят два идентифицируемых параметра 71 и К (в этом уравнении подставлены все промежуточные коэффициенты):
<р(Ь) = К- С-Т.к{^1п(пк4\/Т)/[(пк/Т)2(1 + (лк7У7)2)]х х[ео8(лк-//Т) / (пк/Т) - Т 8ш(лк-//7)] / (пк/Т)} = ¿0,
где С = 4ЩТ-Ь).
Поскольку идентифицируются два параметра, то следует выбрать еще одно граничное условие, чтобы иметь систему двух уравнений для их определения. Таким условием может быть время /3 возвращения судна на исходный курс. Это значит, что в этот момент времени изменение курса обращается в нуль: <^3) = 0. При этом, с целью упрощения, будем считать, что исходный курс равен нулю. Этот момент времени фиксируется в процессе испытаний без всякого труда. Введенное дополнительное условие приводит к следующему соотношению:
К- 2к Ьк [соъ(пЫъ/Т) / (пк/Т) + Т • эш(пк-(3/Т)] / [1 + (пкТ^Т)2] = 0. (10)
Оно является трансцендентным уравнением только для одного из идентифицируемых параметров - 71. Второй параметр К входит в него мультипликативно и просто сокращается. Естественно, что решать его относительно параметра 71 можно только численно, для чего уравнение необходимо преобразовать, выразив явно одно из входящих в него 71. Для этого из суммы выделим первое слагаемое, из которого выразим явно 71 в виде:
71 = -[(Ек*)-(л/7)2(1 + (пТ\/Т)2) / яЩп^/Т) + со$(пГ3/Т)/(п/Т)] / эш^э/Т), (11)
где под символом (2к*) кроется сумма из тех же членов, что и в выражении (10), но суммирование начинается с индекса к, равного 3. Такая форма исходного трансцендентного уравнения допускает простое численное решение методом последовательных приближений Ньютона. После определения параметра 71 легко найти параметр К из уравнения (9) простым делением ¿0 на множитель при К в левой части уравнения.
4. Численная обработка модельных данных
Проведем модельные расчеты по предложенным формулам для характерных случаев маневрирования зигзагом. Как и при расчетах, приведенных в табл. 1, возьмем Т = 75 с, отклонение руля на каждый борт ¿0= 10° и время перевода руля в это положение ^ = 10 с. Это соответствует параметрам маневрирования зигзагом судна среднего тоннажа.
У нас остался в запасе параметр ^ - время возвращения судна на первоначальный курс. Придавая ему различные значения, мы будем получать соответствующие значения идентифицируемых параметров Т\ и К. Пределы возможного изменения достаточно жестки. С одной стороны, (3 не может
быть больше чем (T -t1), т.к. в этот момент угол отклонения от первоначального курса должен стать равным ¿0. С другой стороны, не может быть меньше t1, т.к. в этот момент угол отклонения курса от первоначального значения равен -£0. Программа на любом языке программирования для нахождения Т1 пишется элементарно. Приведем фрагмент такой программы на языке VB5, в котором определяется параметрТ1 на основании уравнения (11) методом последовательных приближений. Фрагмент оформлен в виде функции count_t1(), которая возвращает через специальную пользовательскую структуру ExpTimes искомое значение параметра Т1 с заданной степенью приближения. В функцию передаются следующие аргументы: уже упомянутая структура, поля которой связаны с временными данными -моментами t1, t3, периодом T и значением идентифицируемого параметра, заданная погрешность его определения Eps0, номер последнего учитываемого слагаемого Nobr в сумме (Т,к*) выражения (11); кроме того, сама функция возвращает количество итераций Numblter, реализованных внутри функции, для достижения заданной точности результата.
'определение функции
Dim tt0,eps0 as single Public function count t1(ET as ExpTimes,
Dim nobr, Numblter as integer Eps0 as single, Nobr as integer ) as integer
Public Const pi=3.141592 Dim sums, pikt, eps, Y, tt0, ttl, nmb
'определяем пользовательскую структуру nmb=0
Def type ExpTimes tt0=ET.tt0 ' начальное значение параметра
tl as single DO WHILE eps > Eps0
t3 as single sums = 0
T as single with EN
tt0 as single FOR k = 3 TO Nobr
End type IF k MOD 2 = 0 THEN GOTO nxtl
Dim ET as ExpTimes pikt = pi * k/.T
With ET Y = SIN(pikt*.t1)/pikt~2/(1 + (pikt*tt0)~2)
.tt0 = 0 Y = Y*(COS(pikt*.t3)/pikt + tt0*SIN(pikt*.t3))
.t1=10 sums = sums + Y
.t3=50 nxt1:
.T=7 5 NEXT k
end with nmb = nmb + 1
eps0 = 0.001nobr = 9 Y=(pi/.T)~2*(1+ (pi*tt0/.T)~2)/SIN(pi*.t1/.T)
Y=(Y + .T/pi*COS(pi*.t3/.T))/SIN(pi*.t3/.T) 'обращение к функции count t1 tt1 = (-sums) * Y
NumbIter = count_t1(ET, eps, nobr) eps = ABS(tt1 - tt0)
END tt0 = tt1
LOOP
ET.tt0=tt0 ' возвращаем значение параметра
count t1=nmb ' возвращаем число итераций
End function
Для расчетов были выбраны следующие исходные данные: полупериод процесса Т = 75 с, что соответствует судну среднего тоннажа, 10° - амплитуда отклонений руля ¿0, момент времени t3 изменялся от 10 с до 65 с. Результаты расчетов приводятся в табл. 2 по колонкам: момент, параметры модели Ti и K, отношение K/T1, называемое "послушливостью" рулю, амплитуда отклонения курса от первоначального значения (pnaa, момент времени tmax, когда это значение достигается. Последний есть, фактически, момент времени, в который угловая скорость поворота равна нулю - t4.
Процесс решения уравнения (11) итерационный, потому в представленном отрезке программного кода фиксируется количество итераций для достижения указанной точности в 0.001 с. Процесс этот очень быстро сходится и нигде не понадобилось более 4-х итераций.
Таблица 2
t3, c Tb c K 100 K/T1 <Pmax , ГраД tmax, С
10 -49.30 0.0973 -0.197 -12.7 0
15 -30.50 0.0596 -0.195 -12.7 0
19 -21.91 0.0463 -0.211 -12.7 0
25 -13.14 0.0368 -0.280 -12.7 0
30 -7.55 0.0336 -0.445 -12.7 0
35 -2.50 0.0341 -1.364 -12.7 0
37 -0.50 0.0358 -7.160 -12.7 0
38 0.50 0.0370 7.400 -12.7 0
40 2.50 0.0400 1.600 -12.9 3
45 7.55 0.0501 0.664 -15.3 7
50 1 3.1 4 0.0687 0.523 -19.2 11
55 20.20 0.1117 0.553 -27.7 15
60 30.50 0.2638 0.865 -57.7 20
Хорошо видно, что устойчивость судна в начальный момент маневрирования начинается с значения ^ = 37, 38 с, т.к. она определяется условием Тх > 0. Поэтому все вычисленные значения, соответствующие меньшим величинам /3, мы обсуждать не будем - в этом случае линейный анализ управляемости просто не имеет смысла. При росте /3 растет и величина 71, т.е. увеличивается устойчивость судна на курсе. Руль уже положен на другой борт, а отклонение от первоначального курса все еще растет, именно поэтому 1тах > ^ = 10 (6 колонка табл. 2). Растет при этом и сама амплитуда отклонения курса (5 колонка). В классической работе по управляемости судна (Соболев, 1976) предложена простая формула, связывающая 71 и ?тах в нашем обозначении: 71 « 1.45 ?тах. Наши результаты в принципе похожи на это соотношение, но полного совпадения нет по той причине, что формула Г.В. Соболева получена в предположении мгновенной перекладки руля (^ = 0), а у нас она занимает 10 с. Один из результатов расчета (при 71 = 7.55 с) был представлен выше на рис. 1.
В табл. 2 представлены результаты расчета параметров модели для конкретного значения момента времени ^ = 10 с и величины полупериода Т = 75 с. Для других значений этих параметров весь расчет придется повторить. Но можно указать другой способ, когда мы будем оперировать не абсолютными значениями времени, а их относительными величинами, отнеся их к значению полупериода Т. Тогда в расчетах вместо Т будет фигурировать всегда 1, а вместо ^ - относительная величина Т\ = ^/Т. Так как значение ^ нормировано скоростью поворота руля примерно 1 град/с, то мы в расчетах возьмем такой ряд относительных значений: от 0.05 до 0.20 с шагом в 0.01. При этом само ^ будем считать в соответствии с нормативом равным 10 с. Это означает, что предложенный ряд относительного значения Т\ соответствует следующему ряду значений полупериода Т, с: от 200 до 50 с с переменным шагом. Но самое замечательное при этом, что относительное значение параметра модели 3= 71/Т оказывается зависящим не от гь а только от г3 = /3/Т. Покажем вычисленные значения в табл. 3.
Таблица 3
Г3 = t3/T 0.60 0.65 0.70 0.75 0.80 0.85 0.90
3= Ту/Т 0.1005 0.1547 0.2173 0.2953 0.4011 0.5634 0.9665
Можно пользоваться для обработки данных маневрирования табл. 3, но можно указать простую аппроксимацию этих результатов в виде полиномов 2-ой или 3-ей степеней. Они получены методом наименьших квадратов с помощью пакета MathCad7 и итоги аппроксимации таковы:
3= 1.632 - 5.598г3 + 5.095(г3)2 с СКП 9-10-3 (14)
3= -6.018 + 26.579г3 - 39.646(г3)2 + 20.57(т3)3 c СКП 2-10-3. ( )
С помощью табл. 3 или одного из уравнений (14) параметр модели Номото 71 находится следующим образом. В ходе испытаний фиксируется полупериод Т и момент времени t3, когда отклонение курса становится равным нулю. По их отношению г3 из табл. 3 с интерполяцией или по одной из формул (14) найдем относительное значение параметра модели Л. Умножая его на полупериод Т, получаем значение параметра 71 модели. Для нахождения второго параметра К придется все же использовать формулу (9).
5. Возможные расширения сферы применимости
Обратим внимание на еще один наблюдаемый параметр при выполнении маневра - момент, когда угловая скорость поворота обращается в ноль. Это возможно сделать, если используется самописец курса судна или аксиометр. Назовем этот момент времени t4 и напомним, что в табл. 2 он не наблюдался, а вычислялся. Возможно использовать результат такого наблюдения и найти параметр модели 71 из условия a(t4) = 0. Оно приведет к трансцендентному уравнению, подобному (11):
Т = [(2к**)-(л/Т)2(1 + (пТ1/Т)2) / sin^/T) + sin^/T)] / [(л/Г) cos^/T)], (12)
где означает сумму, которая участвует в выражении (6) для угловой скорости, но суммирование начинается с к = 3. Такое выражение для 7\ даст более медленную сходимость итераций, т.к. угловая скорость есть производная по времени от угла отклонения. За счет дифференцирования каждый к-ый член ряда возрастает в (лк/T) раз, что и замедляет сходимость. Тем не менее, были проведены расчеты с использованием уравнения (12). Теперь число итераций возросло до 7-8, что не играет особой роли при современном быстродействии ВТ. При этом результаты идентификации параметров модели практически не отличаются от результатов с использованием (11).
Но зачем говорить о приеме, который по сходимости явно уступает? Дело в том, что это удобное как по способу наблюдения, так и по вычислительной процедуре краевое условие можно использовать
для решения принципиально новой задачи. А именно, идентифицировать параметры в модели управляемости второго порядка, которую обычно записывают в виде:
T1 T2 d2®/ dt2 + (T1 + T2) da/dt + ® = K S. (13)
В этой модели уже три параметра 7Ь Т2 и К. Следовательно, для их определения необходимо три граничных условия, которыми будут теперь <p(t3,) = 0, <a(t4) = 0 и <р(Т - t1) = ¿>0. Это обобщение метода на модель второго порядка предполагается реализовать в дальнейших работах по проблеме параметрической идентификации обобщенной модели Номото.
6. Заключение
Изложенный в работе новый прием обработки результатов натурного маневрирования типа "зигзаг" с использованием разложения закона перекладки руля в ряд Фурье дал возможность обработать экспериментальные данные по минимальному и легко измеряемому набору характерных времен процесса: времени нулевого отклонения курса или нулевой угловой скорости поворота. Этого достаточно для полной идентификации параметров модели Номото. Более того, получено полуэмпирическое соотношение между параметром модели 71 и моментом времени t3 в относительных величинах. Способ оправдал возлагаемые на него ожидания по упрощению процесса идентификации и показал принципиальную возможность применить тот же прием для идентификации обобщенной модели с большим числом идентифицируемых параметров.
Литература
Nomoto K. and Norrbin N. A review of methods of defining and measuring the manoeuvrability of ships.
ITTC, Manoeuvrability Committee Report, 1969. Nomoto K., Taguchi T. and Hirano S. On the steering qualities of ship. International Shipbuilding Progress, v.4, N 35, p.56-64, 1957.
Асиновский B.A. Использование численного метода для анализа результатов натурных испытаний управляемости судов при неустановившемся движении. Сборник статей молодых научных работников, часть VIII, 1969. Асиновский В.А., Гофман А.Д. Об оценке управляемости судов. Научно-техническое общество
им. акад. А.Н. Крылова. Л., Судостроение, вып.90, 1967. Басин А.М. Ходкость и управляемость судов. М., Транспорт, с.455, 1977.
Кацман Ф.М., Музыкантов Г.М., Шмелев A.B. Эксплуатационные испытания морских судов. М., Транспорт, с.272, 1970.
Соболев Г.В. Управляемость корабля и автоматизация судовождения. Л., Судостроение, с.477, 1976.
