4. Движение цилиндра в жидкости: минимизация работы сил сопротивления

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

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

Положение системы описывается обобщенными координатами $\xi$, $\varphi$. Пусть $v$ - линейная скорость точки захвата, а $\omega$ - угловая скорость цилиндра.

Движение крана осуществляется при помощи управляющей силы $F$, линия действия которой описывается уравнением $y=y_T$. Из всех остальных действующих на кран сил ниже учитывается только сила тяжести. На цилиндрический контейнер действует управляющий момент $U$, сила тяжести, сила Архимеда и гидродинамическая сила. Работа, совершаемая управляющими воздействиями $F$ и $U$, описывается дифференциальным уравнением

$$\dot W=(vF+\omega U)\;.$$

Теперь может быть поставлена следующая задача.

Задача 4.1. Требуется найти управляющие воздействия $F^0(t)$ и
$U^0(t),$ $0\leqslant t\leqslant t_p,$ переводящие с наименьшими затратами энергии $W(t_p)$ цилиндр за заданное время $t_p$ на заданное расстояние $\xi_p$. В начальный и конечный моменты цилиндр имеет вертикальную ориентацию.

Ниже в краткой форме описано решение этой задачи с учетом только энергетических затрат на преодоление силы лобового сопротивления. Полное изложение результатов дано в статье [28].

Рассматриваемая задача является сингулярной и не может быть непосредственно решена при помощи классического вариационного исчисления. Ниже излагается редукция задачи 4.1 к вспомогательной задаче. Последняя не является сингулярной, однако принадлежит к числу задач негладкой оптимизации [26] и требует для своего решения специального подхода. Редукция опирается на то, что, во-первых, движение цилиндра происходит в потенциальном поле силы тяжести. Во-вторых, часть работы расходуется на изменение кинетической энергии цилиндра. Поэтому, поскольку граничные условия заданы, то варируемая часть работы, зависящая от выбора управляющих воздействий, - это энергетические затраты на преодоление силы лобового сопротивления.

В рамках гидродинамических ограничений L1-L3 и гипотезы 1 на движущийся в жидкости цилиндр действуют сила лобового сопротивления и подъемная сила (2.9)-(2.11), где $S$ - площадь проекции цилиндра на плоскость, перпендикулярную вектору его скорости ${\bf V}$. Площадь $S$ однозначно характеризуется углом атаки цилиндра $\alpha$, т.е. наименьшим по абсолютной величине углом между вектором ${\bf V}$ и осью цилиндра (при этом направление отсчета по часовой стрелки считается положительным, поскольку движение цилиндра описывается в левой системе координат). Расчет площади $S$ приводит к формуле ($b=2ld$, $a=\pi d^2/4$)

$$S= V^{-1}\,(b\vert p\vert+a\vert q\vert)\;,$$ (4.1)

где приняты обозначения

$$p=v\cos\varphi-l\omega\;,\quad q=v\sin\varphi\;.$$ (4.2)

С учетом (5.1) выражение для мощности имеет вид

$$\dot A = C_D\,\rho\,{V^2 \over 2}\,(b\vert p\vert+a\vert q\vert)\;, \quad V^2=v^2-2vl\omega\cos\varphi+(l\omega)^2.$$ (4.3)

Направление оси цилиндра определяется вектором

$${\bf E}=-{\bf i}\sin\varphi+{\bf j}\cos\varphi \,.$$ (4.4)

Из определений операции ${\bf a}\to {\bf a}^\bot$ (см. раздел 1), выражений для величин $p$ и $q$ и из формулы (5.4) вытекают соотношения

$$({\bf V},{\bf E})=-q\;,\quad ({\bf V},{\bf E}^\bot)=-p\;,\quad V^2=p^2+q^2\;,$$

которые позволяют записать выражение (2.12) для угла атаки в виде

$$\alpha=-\arctg {p\over q}\;.$$ (4.5)

Представления (5.3) и (5.5) позволяют перейти к следующей вспомогательной задаче динамической оптимизации.

Задача 4.2. Требуется найти функции $v(t)$, $\omega (t)$, минимизирующие терминальный функционал $A(t_k)$ при динамических ограничениях

$$\begin{array}{l} \dot A = C_D(r,-\alpha)\,{\rho V^2\over 2}... ...)=0\;,\quad\dot\varphi =\omega\;,\quad\varphi (0)=0 \end{array}$$ (4.6)

и граничных условиях

$$\xi (t_{\,k\,})=\xi _{\,k\,}\,,\quad \varphi (t_{\,k\,})=0.$$ (4.7)

Угол атаки подсчитывается по формуле (5.5).

Задача имеет две особенности. Во-первых, оптимальный процесс обладает свойством симметрии относительно момента времени $t_k/2$. Поэтому достаточно найти экстремальные управления для некоторых граничных условий $\xi(t_k/2)=\hat{\xi}_k/2$ при свободной ориентации цилиндра в момент времени $t_k/2.$ Во-вторых, в силу отсутствия ограничений на величину управлений задача относится к числу задач классического вариационного исчисления. Но формальной попытке составить уравнения Эйлера-Лагранжа препятствует то, что мощность, а отсюда и гамильтониан не дифференцируемы в ситуации, когда или $p=0$, или $q=0$. Показано, что экстремали находятся среди управляемых траекторий, для которых $q=0$, если $0\leqslant t\leqslant t_1$, $p=0$, если $t_2\leqslant t\leqslant t_3$ и $\omega=0$, если $t\geqslant t_3$. Здесь $t_1\leqslant t_2\leqslant t_3$.

Теорема 4.1. На оптимальном перемещении цилиндра мощность не изменяет своего значения. Производная мощности по скорости точки захвата вне интервалов движения, во время которых $p=v\cos\varphi-l\omega=0$, также сохраняет постоянное значение.

Затем установлено, что существуют экстремали только трех типов. Первая из них соответствует движению цилиндра с сохранением вертикальной ориентации и с постоянной скоростью точки захвата. Уравнения движения в соответствии со второй экстремалью в промежутке времени $[t_1,t_2]$ приведены к виду, удобному для организации вычислительного процесса для нахождения текущих оптимальных значений скоростей $v$ и $\omega$, а именно:

$$\dot\xi = V \frac{\sin\gamma }{\sin\varphi}\,,\quad \dot\varphi = l^{-1} V \frac{\sin(\gamma-\varphi)}{\sin\varphi}\;,$$ (4.8)

$$V = C^{1/3}_2\,C^{-1/3}_D\,N^{-1/3}\,T^{-1/3}\,(\rho/2)^{-1/3}\;,$$ (4.9)

$$\gamma=s\varphi_m+\sigma\arccos T\;,\quad s=s_p\, s_q\;,\qua... ...rm sgn}\nolimits p\;,\quad s_q=\mathop{\rm sgn}\nolimits q\;,$$ (4.10)

$$\sigma(\gamma)=\left\{ \begin{array}{lll} -1\;, & 0\leqsla... ...le {\pi\over 2}<\gamma\leqslant 2\pi\; & \end{array}\right. ,$$ (4.11)

$$\begin{array}{l} \!\!\!\!\!\!\!\!\!\!(2T^2\!+\!1)\cos(\varph... ...hi-s\varphi_m)- T\sin(\varphi-s\varphi_m)\bigr)T\;, \end{array}$$ (4.12)

$$R=\cos^{1/3}\varphi_m\, \left(\frac{\displaystyle C_D(r,{\pi\over 2})}{\displaystyle C_D(r,-\alpha)}\right)^{1/3}.$$ (4.13)

Здесь угол атаки $\alpha=\gamma-\pi/2$, $\varphi_m$ определяется соотношениями $b=N\cos\varphi_m$, $a=N\sin\varphi_m$.

Условия сопряжения экстремальных дуг дают возможность выписать выражения для двух углов $\varphi_{p1}$ и $\varphi_{p2}$, позволяющие определить угол, начиная с которого цилиндр должен двигаться с ориентацией, обеспечивающей равенство $p=0$. Такой режим назван скольжением, так как в ходе него цилиндр движется с нулевым углом атаки. Скольжение описывается системой (5.6) при условии, что мощность постоянна и $p=0$. Такая система интегрируется до конца, при этом

$$\varphi(t)=\arccos e^{- V^0 t/l}\;,\quad \xi(t)= V^0\,t+l\,\ln(1+\sqrt{1-e^{-2 V^0 t/l}}),$$

где $V^0$ - величина скорости цилиндра в скольжении. Если скольжение наступает с момента $t=0$, то мы имеем третью экстремаль.

Если относительное удлинение цилиндра $r=2l/d<10$, то можно выбрать интервал чисел Рейнольдса $10^4<{\rm Re}<5\cdot 10^5$, на котором известны данные по значениям коэффициентов лобового сопротивления цилиндров с образующими как перпендикулярными, так и параллельными потоку (см., например, [12]. Практически, эти значения не зависят от числа Рейнольдса. Эта информация положена в основу интерполирования коэффициентов лобового сопротивления цилиндра по двум переменным - относительному удлинению и углу атаки. Оптимальные скорости $v(t)$ и $w(t)$ ищутся среди тех, которые обеспечивают удержание чисел Рейнольдса на указанном выше интервале. Тогда интерполированные коэффициенты лобового сопротивления - медленно изменяющиеся функции числа Рейнольдса и, следовательно, выполняется гипотеза 2 раздела 1.

Предварительный теоретический анализ позволяет установить интервалы знакопостоянства параметров $s_p$, $s_q$ и $\sigma$. Оказывается, что $s_p=s_q=1$, $\sigma=-1$ для $\varphi<\varphi_p$ и $\sigma=1$ для $\varphi>\varphi_p$, где $\varphi_p=-\varphi_m+\arccos\cos^{1/3}\varphi_m$. Кроме этого, установлено существование двух универсальных констант $\varphi_0$ и $\varphi_1$, которые позволили классифицировать и описать структуру второй экстремали. Если $\varphi_m<\varphi_0$, то сход в режим скольжения происходит по достижению угла $\varphi_{p2}>\varphi_p$. В противном случае это происходит на углу $\varphi_{p1}$. Если $\varphi_m<\varphi_1$, то угол $\varphi_{p1}<\varphi_p$, в противном случае угол $\varphi_{p1}>\varphi_p$. Полученная информация о структуре системы (5.8)-(5.13), являясь полезной, не снимает трудность ее численного интегрирования. Дело в том, что хотя система(5.8)-(5.13) выписана в нормальной форме Коши, она определяет $\dot\xi$, $\dot\varphi$ лишь неявным образом. Это затруднение преодолевается введением для угла $\gamma$ запаздывания по времени, и вычислительный процесс осуществляется по следующей схеме. Пусть в момент $t=t_i$ заданы углы $\varphi(t_i)$ и $\gamma(t_{i-1})$. Тогда можно найти угол атаки $\alpha(t_i)=\gamma(t_{i-1})-\pi/2$, коэффициент лобового сопротивления $C_D(r,-\alpha(t_i))$, коэффициент $R$, корень уравнения (5.12). Затем по формулам (5.10), (5.11) можно подсчитать угол $\gamma(t_i)$, а по формуле (5.9) - скорость перемещения центра масс $V(t_i)$. Эта информация позволит определить правые части системы (5.8) в момент $t_i$. После этого обобщенные координаты цилиндра $\xi (t_{i+1})$, $\varphi (t_{i+1})$ в следующий момент $t=t_{i+1}$ находятся, например, по формуле метода Рунге - Кутта. На рис. 2 представлены графики изменения обобщенных скоростей, угла атаки и скорости центра инерции цилиндра в функции его угловой координаты.

Рис. 2. Графики изменения обобщенных скоростей, угла атаки, скорости
центра инерции цилиндра и корня основного уравнения.

При интегрировании системы (5.8)-(5.13) необходимо на каждом шаге решать уравнение (5.12). Главная особенность этого уравнения состоит в неединственности его решений. Разработан алгоритм, который обеспечивает удержание нужного корня на нужной ветви решений уравнения (5.12) и в случае исчезновения этой ветви переводит его на другую ветвь. Чтобы найти нужный корень на очередном шаге интегрирования системы (без ``срыва'' на чужую ветвь), применяется специальная процедура. Она состоит в нахождении для искомого корня симметричного интервала. Эти и другие алгоритмические особенности, сопровождающие решение рассматриваемой задачи, детально обсуждаются в [28].

Вычислительный эксперимент позволил сделать вывод, что на всех удлинениях $r$ из интервала (0.4, 4.2) вторая экстремаль не может дать результат лучше, чем третья. Система (5.6) интегрируется в квадратурах как в случае движения цилиндра в соответствии с первой экстремалью, так и с третьей. Это дает возможность сравнения первой и третьей экстремалей уже аналитическими методами.

В результате установлены следующие факты. Если относительное удлинение цилиндра превышает величину $(L/d)_{cr}$ и отношение $\xi_{k}/l$ больше, чем некоторое критическое значение $(\xi_{k}/l)_{cr}$, зависящее от $L/d$, тогда третья экстремаль является оптимальной. В противном случае оптимальным является движение цилиндра с сохранением вертикальной ориентации и с постоянной скоростью. Здесь $(L/d)_{cr}$ и $(\xi_{k}/l)_{cr}$ - корни некоторых трансцендентных уравнений. Значение $(L/d)_{cr}$ как функция относительного удлинения цилиндра неограниченно и монотонно возрастает при $(L/d)\to(L/d)_{cr}+0$.

Исходную задачу можно решить как обратную задачу динамики. Действительно, информация об оптимальных управлениях в задаче 4.2 является по сути дела информацией об оптимальных фазовых траекториях в исходной задаче 4.1. Следовательно, достаточно подставить оптимальные $v(t)$ и $\omega (t)$ в уравнения движения цилиндра:

$$\begin{array}{l} ml\dot v\cos\varphi + (ml^2+J)\dot\omega = ... ...mega\cos\varphi +ml\omega^2\sin\varphi = Q_{\xi}\;, \end{array}$$

где $J$ - момент инерции цилиндра,

$$\begin{array}{l} Q_{\varphi}=-mgl\sin\varphi -l\bigl((D_x+D_... ...hi )\bigl)+M_G+U\;,\\ [7pt] Q_{\xi}=D_x+D_x^l +F\;, \end{array}$$

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