Моделирование движения стаи рыб в однородно стратифицированной жидкости

Задача поиска стай промысловых рыб представляет собой большой практический интерес. Поиск косяков рыбы может осуществляться непосредственно, когда возможно идентифицировать скопление рыбы гидро-акустическими методами по анализу отраженного стаей звукового сигнала [1], в оптическом диапазоне [2] или по косвенным признакам, когда наблюдения ведутся за ключевыми параметрами среды обитания рыбы, температурой и количеством растворенного кислорода [3], а также по более сложным моделям, включающим в себя данные зондирования температуры поверхности, солености, толщины слоя перемешивания, концентрации планктона [4, 5]. Когда скопление рыбы найдено, в эту область направляются рыболовецкие средства. В настоящей статье проводится моделирование движения стаи рыб в толще океана. Вышедшее на поверхность возмущение от движущейся стаи рыб может быть затем зарегистрировано средствами аэрокосмического зондирования.

• 1.    Постановка задачи

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

Рис. 1 . Расчетная область Q. Схематически показана стая из 27 рыб, центр масс стаи находится в точке x cm (t)

Для описания движения жидкости будем использовать приближение Бусси-

неска [11,12], запишем уравнение на внутренний потенциал ^(x, t) в области Q =

[x min , x max

] ^х [У тт , y max

] х [-H,H ]

. d ² ^

А оТГ + N ² ^ xy ^ = f ( x ,t), x ^E Q ∂t

с граничным условием

^ = 0, x G dQ,

где N - частота плавучести. Считается, что область Q (см. рис. 1) достаточно велика, так что граничное условие (2) и отраженная от границы волна возмущения не влияют на результаты расчета вблизи источника возмущения. В начальный момент времени

#=о = 0.

Уравнение (1) является уравнением соболевского типа [13], точные утверждения о разрешимости задачи (1) - (3) для случая N = const можно найти в [14].

Пусть стая состоит из M одинаковых рыб, тогда правая часть (1) представляет собой сумму ^А=1 B5 i (x i (t)), где £( • ) - 5-функция Дирака. Для упрощения численного счета будем задавать f ( x ,t) в виде

м f (x,t) = ^ fi(x - xi(t)), i=1

где X i (t) - положение i-й рыбы в момент времени t, функция f i определяется константами A _i , B _i :

^f i ^(x,y^,z) =

B i A i „-А? ( x ²+ y ² +z ² ) ^П3

Заметим, что f f i (x}d x = B i .

R ³

Скорость жидкости v = (v _x ,v y ,v _z ) может быть найдена как

v = д²^ ^v + N ^{2 v} *y) ^ф,                      (6)

давление p(x, t):

P(x,t} = -Po      + N2^) -d^(x,t),

∂t2

вертикальное смещение жидкости Z( x ,t) :

∂∂

Z(x,t) = oilZ^t)'

Для решения задачи (1) – (3) на языке программирования С++ была написана компьютерная программа, реализующая неявную разностную схему на равномерной расчетной сетке [15], где для решения системы линейных уравнений на каждом шаге по времени использовался обобщенный метод минимальных невязок [16]. Рассчитывались поля скорости (6), давления (7), вертикального смещения (8) жидкости. С использованием этой программы проводилось моделирование движения одиночного источника возмущения в однородной стратифицированной жидкости, при этом точечный массовый источник возмущения двигался равномерно горизонтально [17, 18] или под углом к горизонту [15, 19], корректность расчетов верифицировалась сравнением с известными асимптотическими решениями [20, 21], c результатами экспериментов [22]. Модифицированный вариант программы использовался для моделирования возмущения с условием свободной поверхности на верхней границы области Q; в этом случае решалась задача типа (1) для вертикальной компоненты скорости vz в области Q, а затем задача Коши для обыкновенного дифференциального уравнения в каждой точке вертикальной границы z = H области Q, см. [23]. Вообще говоря, реализованный метод расчета позволяет проводить моделирование для изменяющейся по высоте частоты плавучести, N = N(z), когда частота плавучести изменяется скачком [24] или непрерывно [15], однако далее в настоящей работе будем считать, что жидкость однородно стратифицированная, N = const. Исходный код расчетной программы доступен в интернете по адресу:

1.1.    Численное моделирование движения стаи рыб

В настоящей главе численно исследуется возмущение от движущейся стаи рыб, его характерные отличия от возмущения, создаваемого отдельным точечным массовым источником. В частности, для интересов обнаружения стай промысловых рыб, представляет интерес изучение полей возмущения вблизи стаи рыб, а также возмущения, появившегося в результате резкого маневра стаи. Такие исследования требуют численного решения задачи (1) – (4), они не могут быть проведены аналитически или с использованием известных асимптотических формул [20, 21], верных в дальней зоне для устоявшихся присоединенных волн.

Будем считать, что стая состоит из M одинаковых рыб, расположенных в узлах правильной кубической сетки с шагом р. Схема такой стаи для M = 27 показана на рис. 1. Функция f (x, t) задающая правую часть (1) определяется формулами (4), (5). Далее в численных расчетах жидкость считалась равномерно стратифицированной, частота плавучести N =1; значение констант в (5) принималось равным A i = B i = 1, i = 1,...,M .

1.1.1.    Движение с постоянной горизонтальной скоростью

В случае, когда все особи движутся с одинаковой и постоянной по величине и направлению скоростью U = (a,b, c), для положения i-й рыбы в момент времени t имеем:

X i (t) = x 0 + Ut,                                    (9)

где x ⁱ ₀ – начальное положение i -й рыбы.

z 0

-40

80     120    160

0.18

0.16

0.14

0.12

0.1

0.08

0.06

0.04

0.02 0

Рис. 2 . Функция f , задающая правую часть (1), в вертикальном сечении у = 0 (a), нормированное вертикальное смещение жидкости Z (x, 0,z) в вертикальном сечении у = 0 (b), горизонтальное сечение функции вертикального смещения жидкости

Z (x, у, h) , сечение находятся на расстоянии h = 20m от плоскости движения центра масс стаи (c). От момента начала движения прошло t = 82, 6 с

Рассмотрим пример численного моделирования горизонтального прямолинейного движения. Пусть стая из M = 27 рыб движется горизонтально в положительном направлении оси Ox со скоростью U = (1, 0, 0), см. рис. 1. Расстояние между любыми двумя ближайшими рыбами p = 5 м, то есть стая имеет размер 10 м х 10 м х 10 м. В начальный момент времени центр масс стаи x ^cm = M ^ ^1 x i , находится в точке с координатами ( - 40, 0, 0). Схема стаи и функция f (x, t) показаны на рис. 1 и рис. 2 а соответственно, результаты численного решения задачи (1) – (5), (9) приведены на рис. 2 б и рис. 2 в. На рис. 2 б показано нормированное вертикальное смещение жидкости Z(x,y,z,t) = N2nB 1 Z(x,y,z,t) в вертикальном сечении у = 0, проходящем через центр стаи. Видна картина обтекания рыб стаи, находящихся на уровнях глубины z = — p, z = 0, z = p. На рис. 2 в показано вертикальное смещение жидкости Z (x, у, 20) в горизонтальной плоскости, находящейся на уровне z = 20, то есть на расстоянии 15 м от стаи по вертикали. В картине возмущения на рис. 2 в видно влияние формы стаи по сравнению с картиной возмущения от одиночного точечного источника [20]. Один такой расчет, то есть моделирование движения стаи из точки ( — 40, 0, 0) до момента ее выхода на границу области Q, занимал примерно трое суток на одном вычислительном узле (два процессора Intel Xeon 3 ГГц, 192 Мб ОЗУ) суперкомпьютера МВС-10П Центра коллективного пользования вычислительными ресурсами МСЦ НИЦ Курчатовский институт. Для расчета использовалась разностная сетка размером 300 х 300 х 300 точек, было сделано 90 шагов по времени.

1.1.2.    Движение в горизонтальной плоскости с поворотом на 90°

Пусть стая рыб начинает движение в момент времени t = 0 так, что каждая рыба движется в положительном направлении оси Ox со скоростью U 1 = (a1, 0, 0). Центр масс стаи в начальный момент времени находится в точке с координатами x 0m = (x 0 ^m , 0, 0), x 0 ^m < 0. В момент времени, когда центр масс достигает плоскости x = 0, скорость каждой рыбы мгновенно становится равной U ₂ = (0, b ₂ , 0). Таким образом, скорость изменяется со временем следующим образом:

ТТМ   Г U1 = ^{( a} 1 , ⁰ , ⁰⁾ , ⁰ <  t <  ^{— x} 0m ^/a 1 ,

^{U ( t )} = t U 2 = (0,b ₂ , 0), t> —x0 ^m /ai .

Далее показаны результаты моделирования поворота на 90° для скоростей a1 = b2 = 1 м/с. Стая состоит из M = 27 одинаковых рыб, расстояние между рыбами p = 5 м. В начальный момент времени центр масс стаи x0m = x04 находится в точке с координатами (—40 м, 0, 0). Схема движения показана на рис. 3. Результаты численного решения задачи (1) – (5), (9), (10) показаны на рис. 4. Показано горизонтальное сечение z =1 функции вертикального смещения жидкости Z(х,У, 1) в момент t = 39, 8 c непосредственно перед поворотом (рис. 4 а), в момент времени t = 63, 7 c (рис. 4 б), в момент времени t = 132, 3 c (рис. 4 в). Из рисунков видно, что сначала за движущейся стаей формируется картина присоединенных волн. После скачкообразного изменения направления скорости происходит деформация картины возмущения в горизонтальной плоскости, а затем, со временем, за продолжающей движение вдоль оси Oy стаей формируется классическая картина присоединенных волн; одновременно с этим продолжается (с меньшей интенсивностью) распространение возмущения по первоначальному направлению, а именно вдоль оси Ox.

U 2

U 1

U 2

—►

x

Рис. 3. Схема движения стаи из M = 27 одинаковых рыб. Движение происходит в горизонтальной плоскости z = 0. Сначала стая движется вдоль оси Ox. В момент, когда центр стаи оказывается в начале координат, направление скорости изменяется на 90 градусов, далее стая движется вдоль оси Oy у о

-40

-80

0.08

0.06

0.04

0.02

-0.02

-0.04

-0.06

-0.08

-0.1

-80 -40   0   40  80

x

у 0

-40

-80

-80

-40   0

x

0.08

0.06

0.04

0.02

-0.02

-0.04

-0.06

-0.08

-0.1

Рис. 4 . Горизонтальное сечение z = 1 функции вертикального смещения жидкости Z (x, у, 1) в момент непосредственно перед поворотом t = 39, 8 с, в момент времени t = 63, 7 с, в момент времени t = 132, 3 с (слева направо)

1.1.3.    Движение рыб по синусоидам

В расчетах, результаты которых приведены выше стая рыб двигалась поступательно как целое, с постоянным расстоянием между рыбами, причем каждая отдельная рыба двигалась прямолинейно с постоянной скоростью. Однако в природе каждая рыба движется по изогнутой траектории, а расстояние между рыбами постоянно меняется. Промоделируем такое более сложное движение. Пусть координаты i -й рыбы изменяются по следующему закону, i = 1,..., M :

xi (t) = x0i + at, yi(t) = yoi + dsin (2f t + фi^ ,                            (11)

Z i (t) = Zo i .

то есть, траектория i-й рыбы лежит в горизонтальной плоскости z i = zo i и представляет собой сдвинутую на случайную величину ф _i Е [0, 2п] синусоиду с периодом l f и амплитудой d, при этом проекция скорости рыбы на ось Ox постоянна и равна a.

На рис. 5 показаны результаты расчета движения стаи рыб из M = 125 особей. Каждая рыба движется по закону (11) с параметрами l f = 1 м, d = 0, 2 м, центр масс стаи движется в положительном направлении оси Ox со средней скоростью a = 1 м/с. Точки (x gi , y oi , z gi ), i = 1,..., 125 (то есть координаты рыб в начальный момент времени, если бы все случайные сдвиги фаз φ _i были равны нулю) находятся в вершинах кубической сетки 5 х 5 х 5 с шагом 1 м и с центром в точке с координатами ( - 40, 0, 0).

z

-30

-4

-2

-6

x

-30

-8

z

-30

0.8

0.6

0.4

0.2

-0.2

-30      0      30

У 0

-30

0.06

0.04

0.02

-0.02

-0.04

-0.06

-30      0      30

x

Рис. 5 . Моделирование синусоидального движения рыб в стае. Показаны нормированное вертикальное смещение Z (x, 0, z) в вертикальном сечении y = 0 м (a), горизонтальная компонента скорости v _y в вертикальной плоскости z = 1 м (b) и в горизонтальной плоскости z = 10 м (c) в момент времени t = 28, 8 с после начала движения

Сравнивались результаты расчетов для случаев движения рыб по синусоидам (11) и их прямолинейного движения, то есть когда d = 0 в (11). Сравнение результатов проводилось в горизонтальном сечении z = z ₀ , —H < z ₀ < H , в области Q _z 0 = [x min ,x m_ax ] х [y min ,y _max ]. В качестве меры отличия использовалась величина e(g s ,g ^r ):

^ ⁽ g ^{s ,}g ^r ) =

|| g ^{s ( x )} — g ^{r ( x ) ||} L 2 (Q z0 ) || g r ^{( x ) ||} L 2 (Q z ₀ )

где || g || _L 2 (Q _z 0 ) = у Q _o g(x, У, z ₀ ) ² dxdy, g s — некоторая гидромеханическая характеристика для синусоидального движения рыб, g ^r – та же характеристика, но для случая прямолинейного движения рыб. Расчеты показали, что для горизонтальных компонент скорости, для вертикального смещения жидкости величины отличия e(v X ,v X ), e(v y ,v y ), e(Z X ,Z X ) составляют не более 1%. Таким образом, при рассмотренных нами значениях параметров результирующие поля возмущений для случаев синусоидального и прямолинейного движения рыб в стае отличаются несущественно.

1.1.4.    Исследование зависимости характера возмущения от расстояния между рыбами в стае

Изучим влияние расстояния между рыбами на картину возмущения при горизонтальном поступательном движении стаи. Пусть стая рыб состоит из M = 125 одинаковых рыб, находящихся в начальный момент времени в точках с координатами xiJ,k = (хГ + i • p,yr + j • p, ^Г + k • p), i,j,k = —2, —1,0,1,2, где (x0m, yCm,zCm) — координаты центра масс стаи в момент времени t = 0, а именно, xcom = —40, yocm = zocm = 0. Каждая рыба движется с горизонтальной скоростью

U _i = (1, 0, 0), i = 1, . . . , 125. Результаты моделирования для различных значений расстояния между рыбами p показаны на рис. 6. Показаны возмущение поля горизонтальной компоненты скорости v x в горизонтальных сечениях z = 10 м (рис. 6 ж, 6 з, 6 и), z = 20м (рис. 6 г, 6 д, 6 е) и z = 70м (рис. 6 а, 6 б, 6 в) в момент времени t = 43, 8 с. Расстояние между рыбами p принимало значения 1 м (рис. 6 ж, 6 г, 6 а), 1,5 м (рис. 6 з, 6 д, 6 б) и 2 м (рис. 6 и, 6 е, 6 в). Результаты расчетов показывают, что вблизи движущейся стаи характер возмущения гидродинамических полей существенно зависит от расстояния между рыбами, тогда как на достаточно большом расстоянии от косяка рыб картина возмущения уже не зависит от расстояния между рыбами и близка к картине устоявшихся присоединенных волн в дальней зоне [20].

y 0

-40

-80

-80

-40

0.02

0.015

x y 0

-40

-80

0.01

0.005

-0.005

-80

-40

x y0

-40

-40

-80

x

0.06

0.04

x

У 0

-40

-80

x

-40

-80

-0.01

-0.015

-0.02

-80

Рис. 6 . Возмущение поля горизонтальной компоненты скорости v x в плоскостях z = 10 м (нижний ряд), z = 20 м (средний ряд) и z = 70 м (верхний ряд) в момент времени t = 43, 8 с. Движение стаи из 125 особей происходит в горизонтальной плоскости. Каждая рыба движется прямолинейно со скоростью 1 м/с вдоль оси Ox. Расстояние между рыбами p равно 1, 1,5, 2 (слева направо)

0.02

-0.02

-0.04

x

-40

-80

-0.06

Заключение

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

Было проведено моделирование движения стаи промысловых рыб из 27 или 125 особей. Стая двигалась горизонтально прямолинейно или с поворотом на 90 градусов. Величина скорости движения каждой рыбы была постоянна. В случае движения с поворотом, в момент поворота направление скорости мгновенно менялось на 90 градусов в горизонтальной плоскости. Было показано, что вблизи стаи и в момент изменения скорости движения (в момент начала движения, в момент поворота) картина возмущения существенно отличается от картины устоявшихся присоединенных волн в дальней зоне от одиночного источника возмущений.

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

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

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

Исследование выполнено за счет гранта Российского научного фонда №24-6100025,