A simple mechanical model of turbulence

Full Text

1. Предварительные сведения
Исследование турбулентности является одной из важнейших задач как прикладной, так и
фундаментальной науки. Обусловливается это широким распространением указанного явления
в природе, особенно в прикладных механических задачах, связанных с течением жидкости и
газа, а также в разнообразных атмосферных процессах и явлениях. Несмотря на стремительное
развитие вычислительной техники и разработку численных методов с соответствующим матема-
тическим обеспечением, моделирование турбулентных течений остается одной из сложнейших
проблем механики жидкости и газа. Традиционно для описания турбулентности используется
три различных подхода — статический, структурный и динамический. Однако каждый из
них не дает общего математического способа описания этого физического явления, так как
области применимости каждого из методов существенно различны [1–4]. Поэтому создание
новых подходов к описанию и моделированию турбулентного движения представляется важным
и востребованным.
Кратко опишем основные подходы к моделированию турбулентных движений. К первому
из них относится метод, основанный на использовании усредненных по Рейнольдсу уравнений
Навье–Стокса (Reynolds Averaged Navier–Stokes — RANS) с замыканием с помощью той или иной
полуэмпирической модели турбулентности. Ко второй категории относятся два классических
вихреразрешающих подхода, а именно прямое численное моделирование турбулентности (Direct
Numerical Simulation — DNS) и моделирование крупных вихрей (Large Eddy Simulation —
LES). Наконец, третья категория включает гибридные подходы, опирающиеся на совместное
использование RANS и LES в различных областях потока [5–9].
Семенов М.Е., Толкачев А.В., Канищева О.И. Механическая модель турбулентности
Semenov M.E., Tolkachev A.V., Kanishcheva O.I. A simple mechanical model of turbulence 64 из 75
Исчерпывающее и точное определение турбулентности дать крайне тяжело. Возможно
только выделить ряд свойств, которые будут описывать это явление. Одним из таких свойств,
характеризующих турбулентный поток, является ”завихренность” (растяжение вихрей). Вихри,
которые образуются в турбулентном потоке, сохраняются достаточно долгое время благодаря
механизму каскадного переноса энергии. Когда большие и неустойчивые завихрения, ”отбирая”
энергию у основного потока, временно сохраняют ее, после чего растягиваются и разбиваются
на несколько меньших, и так до тех пор, пока не рассеют кинетическую энергию в тепло. Таким
образом, осуществляется направленная передача энергии от более крупных вихрей к более
мелким, формирующая их иерархию в потоке [10; 11].
Особого внимания при рассмотрении турбулентного течения заслуживают мелкомасштабные
вихри, структура которых подразумевает их однородность и изотропность. Теория локально-
изотропной и однородной турбулентности представляет собой достаточно разработанный раздел
теории турбулентности. Принципиальные результаты в этой области получены А.Н. Колмого-
ровым [12], который предложил ”закон пяти третей”, играющий на сегодняшний день ключевую
роль при разработке современных подходов к моделированию турбулентных течений и их
верификации.
Другим примером сложного нелинейного явления, играющего значительную роль при проек-
тировании и исследовании сложных динамических систем, является гистерезис. Нелинейности
гистерезисной природы присущи многим физическим, биологическим и экономическим процес-
сам и явлениям. Во многом это либо обусловливается их внутренней структурой (магнитный
гистерезис), либо является следствием их динамических особенностей (рулевой люфт). Строгое
математическое описание гистерезиса было дано М.А. Красносельским и А.В. Покровским в
монографии [13]. В ней гистерезисные нелинейности трактуются как операторы или преоб-
разователи с соответствующими пространствами состояний. Выход такого преобразователя
будет зависеть не только от мгновенного значения входа, но и от состояния преобразователя в
предшествующие моменты времени. Разработанные в [13] модели позволяют использовать их
в уравнениях динамики систем с носителями гистерезисных свойств. При этом модели соот-
ветствующих систем сводятся к операторно-дифференциальным уравнениям, для которых к
настоящему времени доказаны теоремы существования и единственности [14]. Данным задачам
посвящено множество работ, среди которых отметим публикации [15; 16].
В настоящее время для моделирования динамики систем с гистерезисом, как правило, ис-
пользуются два подхода: конструктивный подход (неидеальное реле, преобразователь Прейзаха,
модель Ишлинского [17–22]) и феноменологический подход (модель Боука — Вена, Айвана,
Дьюема [23–27]). Каждый из них зарекомендовал себя в различных областях, например, опера-
тор Прейсаха находит широкое применение при моделировании систем с ферромагнитными
элементами, а также в гидрологических моделях проникновения осадков в почву [28]. Среди
феноменологических моделей наиболее популярным является модель Боука – Вена. Эта модель
формализуется посредством двух соотношений: одного алгебраического и дифференциального
уравнения. Модель Боука – Вена является удобным инструментом для формализации гистере-
зисных зависимостей, особенно в ситуации, когда гистерезисное звено является частью сложной
системы. Применительно к проблемам турбулентных течений отметим, что в традиционных
подходах гистерезисное трение к настоящему времени не рассматривалось, хотя механизмы
диссипации энергии в каскадах турбулентных вихрей допускают ”гистерезисную” трактовку.
Применительно к тематике настоящей работы отметим такое явление, как аэродинамический
гистерезис, его детальное описание можно найти в статье [29].
В настоящей статье исследуется перенос энергии в механической системе, состоящей из свя-
занных осцилляторов, соединенных между собой иерархическим образом. Указанная система
является достаточно простой моделью переноса энергии в турбулентной среде. Каждый осцил-
лятор в системе, кроме первого, связан со своими соседями сильно нелинейными кубическими
пружинами. Помимо пружин в связи между осцилляторами включены вязкие и гистерезисные
демпферы, которые формализуются с помощью феноменологической модели Боука–Вена. В эту
систему гистерезис добавляется как естественный элемент управления и стабилизации. Это
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 63–75
Vestnik of Samara University. Natural Science Series 2024, vol. 3, no. 30, pp. 63–75 65 из 75
исследование помогает найти некоторые общие законы передачи энергии между осцилляторами
и предлагает новый способ изучения турбулентности жидкости без использования системы
Навье–Стокса.
2. Описание системы и ее верификация
В природе и в технических устройствах часто наблюдаются разнообразные явления, вклю-
чающие передачу энергии между элементами в виде каскада-трансляции. Эти элементы при
проектировании находятся по отношению друг к другу на различных уровнях или масштабах,
формируя таким образом иерархию, в которой осуществляется каскад. Такой перенос энергии
описывает передачу преимущественно с более крупных масштабов на более малые, например,
как для однородного изотропного движения в гидродинамике. Часто в качестве модели для
исследования подобного явления используют нелинейные цепочки осцилляторов [30–32].
Цепочка нелинейных осцилляторов является универсальной и относительно простой моде-
лью, способной описать сложные нелинейные процессы. Простейшей из таких моделей является
модель межатомного взаимодействия в кристаллической решетке. Отметим, что в подавляющем
количестве публикаций, посвященных связанным осцилляторам, используются упругие взаимо-
действия. Однако реальная природа межатомных связей зачастую требует учета принципиально
иных форм взаимодействия, а именно взаимодействия гистерезисной природы. Поэтому пред-
ставляется важной задача изучения динамики гистерезисно связанных осцилляторов. В этой
связи отметим близкие к тематике настоящей статьи работы [33; 34].
Еще одно важное приложение динамики связанных осцилляторов обусловливается ис-
пользованием соответствующей модели для описания турбулентных процессов. В настоящей
работе предлагается использовать в качестве простейшей модели турбулентных движений си-
стему обыкновенных дифференциальных уравнений, описывающую в том числе и колебания
связанных осцилляторов. При этом для старого обоснования такой замены требуется прове-
сти верификацию модели, сравнивая полученный экспериментально энергетический спектр с
законом Колмогорова.
Рассмотрим сначала механическую модель исследуемой системы, которая представлена
на рис. 2.1. Она состоит из осцилляторов, каждый из которых кроме первого x0 связан со
своими соседями xn нелинейными пружинами (с потенциалом knx4n
/4), где kn — жесткость
нелинейной пружины для n = 1, 2, . . . , N − 1. Помимо пружин в связи между осцилляторами
включены вязкие и гистерезисные демпферы с коэффициентами l и (1− α)knD соответственно.
Гистерезисный демпфер формализуется с помощью феноменологической модели Боука–Вена,
имеющей простую реализацию в виде обыкновенного дифференциального уравнения.
Рис 2.1. Исследуемая механическая система
Fig 2.1. The considered mechanical system
Семенов М.Е., Толкачев А.В., Канищева О.И. Механическая модель турбулентности
Semenov M.E., Tolkachev A.V., Kanishcheva O.I. A simple mechanical model of turbulence 66 из 75
Уравнения движения для механической системы на рис. 2.1 будут иметь вид:


m0 ¨ x0 + αk0x0 + l0 ˙ x0 + (1 − α)k0Dzl0
+ αk1(x0 − x1)3 + l1(x˙0 − x˙1)+
+(1 − α)k0Dzr0
= f cosΩt,
. . .
mnx¨n + αkn(xn − xn−1)3 + ln(x˙n − x˙n−1) + (1 − α)knDzl
n+
+αkn+1(xn − xn+1)3 + ln+1(x˙1 − x˙n+1) + (1 − α)k0Dzr
n = 0,
. . .
mN−1x¨N−1 + αkN−1(xN−1 − xN−2)3 + lN−1(x˙N−1 − x˙N−2)+
+(1 − α)kN−1Dzl
N−1 = 0,
(2.1)
где zl
n и zr
n — гистерезисные переменные, являющиеся решением следующих дифференциальных
уравнений (2.2) для левой и правой (l и r ) стороны относительно n-го осциллятора.
Имеем
z˙l
n = D−1
A − |zl
n|η(βsign
????
zl
n(x˙n − x˙n−1)

+ γ)

(x˙n − x˙n−1),
z˙r
n = D−1 [A − |zr
n|η(βsign (z˙r
n(x˙n − x˙n+1)) + γ)] (x˙n − x˙n+1). (2.2)
В уравнениях (2.2) α, β, γ, η — безразмерные гистерезисные параметры, которые определяют
форму и размер петли гистерезиса согласно подходу модели Боука–Вена. При исследовании
системы примем следующие значения: A = 1, β = γ = 0.5, η = 4 . Для других же гистерезисных
параметров в системе (2.1) значения будут: D = 1 и α = 0.5.
Параметры mn — масса осциллятора, xn — его смещение относительно положения равновесия,
x˙n — скорость и x¨n — ускорение для соответствующих элементов в цепочке.
Зададим параметры так, чтобы в системе (2.1) была иерархия масштабов в направлении
слева направо. Для этого пусть у первого линейного осциллятора m0 = 1, k0 = 1, l0 = 0.01.
Далее за ним у второго нелинейного осциллятора m1 = 0.01, k1 = 0.1, l0 = 0.001. Параметры
остальных осцилляторов в цепи вплоть до N − 1 задаются по правилам:
mn
mn+1
= 6,
kn
mk+1
= 45,
ln
ll+1
= 2,
где n = 0, 1, 2, ...N − 1.
Теперь проведем верификацию модели. Энергетический спектр ˆE (κ) рассчитаем с помощью
отношения
ˆE
(κ) =
¯En
¯E
,
где ¯En(κ) — усредненная по периоду T энергия колебаний для n-го осциллятора цепи:


¯E
0 = 1T
Rt0
t0+T
???? 1
2m0x˙2
0 + 1
2 k0x2
0 + 18
k1(x0 − x1)4
dt,
. . .
¯E
n = 1T
Rt0
t0+T
???? 1
2mnx˙2n
+ 1
8 kn(xn − xn−1)4 + 18
kn+1(xn − xn+1)4
dt,
. . .
¯E
N−1 = 1T
Rt0
t0+T
???? 12
mN−1x˙2
N−1 + 18
kN−1(xN−1 − xN−2)4
dt,
¯E
— усредненная полная энергия системы для всех осцилляторов, которая рассчитывается
следующим образом:
¯E
=
N−1
Σ
n=0
¯E
n,
а κ = 1
mn — волновое число для n = 0, 1, 2, ...N − 1. В итоге общая математическая схема
верификации модели будет выражаться:
ˆE
(κ) ∼ κ−53
. (2.3)
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 63–75
Vestnik of Samara University. Natural Science Series 2024, vol. 3, no. 30, pp. 63–75 67 из 75
Результат расчета энергетического спектра системы (2.1) приведен на рис. 2.1 в лога-
рифмическом масштабе. Черными сплошными линиями обозначены энергетические спектры,
полученные без учета (при α = 1) гистерезиса в системе, красным обозначены результаты,
полученные с учетом гистерезиса, а синяя линия — это ”закон пяти третей” (2.3).
Моделирование осуществлялось для пятнадцати осцилляторов (N = 15) в условиях резо-
нансного внешнего воздействия на первый элемент в цепи (Ω = 1) при различных амплитудах
( f = 0.001, 0.0025, 0.005 и 0.01). Решение исходной системы (2.1) осуществлялось методом
Рунге–Кутта 4-5-го порядка при нулевых начальных условиях.
Рисунок 2.2 иллюстрирует, как энергия первого осциллятора транслируется остальным ос-
цилляторам цепи и убывает по достаточно сложному закону. Видно, что при высокой амплитуде
воздействия каскад энергии гораздо более интенсивен, а для осцилляторов меньшего масштаба
по наклону он в достаточной степени соответствует гипотезе Колмогорова. Так при ампли-
туде внешнего воздействия кривая спектра (красная линия с квадратом) и соответствующая
логарифмическая прямая ”закона пяти третей” находятся в высокой степени соответствия.
101 103 105 107 109 1011
1/mN−1
10−21
10−18
10−15
10−12
10−9
10−6
10−3
100
EN−1
f=0.001
f=0.0025
f=0.005
f=0.01
f=0.001 with hysteresis
f=0.0025 with hysteresis
f=0.005 with hysteresis
f=0.01 with hysteresis
-5/3 law
Рис 2.2. Энергетический спектр системы (2.1)
Fig 2.2. Energy spectrum of the system (2.1)
3. Особенности динамики исследуемой системы
Изучим подробнее динамику системы. Для этого рассчитаем амплитудно-частотные харак-
теристики у отдельных осцилляторов в цепи с помощью метода ”сканирования” частотой. Суть
данного метода заключается в том, что при компьютерном моделировании системы (2.1) изме-
няется значение частоты Ω внешнего воздействия последовательно, без прерывания процесса
вынужденных колебаний (”сканирование частотой”). Иными словами, изменяя последовательно
частоту Ωi, c каждой i-й итерацией значения для начальных условий при следующей Ωi+1
следует брать не нулевыми, а с того момента времени, когда прервалось моделирование на
текущем i-м шаге.
Особенности динамики без учета гистерезиса подробно изучены в работе, где показано,
что в зависимости от амплитуды внешнего воздействия в цепи наблюдаются различные типы
синхронизации колебаний [35; 36]. Также авторы отмечают, что при высоких значениях f и
частоте Ω = 1 система реагирует хаотично, и значительная часть энергии передается по цепи
в виде энергетического каскада-трансляции. При таких условиях существуют хаотические
области, в которых все нелинейные осцилляторы одновременно возбуждаются до хаотического
Семенов М.Е., Толкачев А.В., Канищева О.И. Механическая модель турбулентности
Semenov M.E., Tolkachev A.V., Kanishcheva O.I. A simple mechanical model of turbulence 68 из 75
состояния (активируются), когда частота Ω превышает нижний критический порог; таким
образом, в системе происходит хаотическая синхронизация.
Результаты расчета амплитудно-частотных характеристик для второго, четвертого и шестого
осциллятора в цепи в условиях гистерезисного демпфирования представлены на рис. 3.1 (a,b,c),
когда на первый осциллятор в цепочке оказывается внешняя сила с амплитудами f = 0.001, 0.005
и 0.01.
Ниже (рис. 3.1, с) проведена индентификация динамических режимов для всех элементов
в цепочке при частоте внешнего воздействия Ω = 0.99. Для этого рассчитывался старший
показатель Ляпунова Λ методом Розенштейна [37] для амплитуд xn каждого осциллятора из
временного окна t ∈ (0, 4000) .
0.6 0.7 0.8 0.9 1.0 1.1 1.2

0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
Amplitude
f=0.001 up
f=0.001 down
f=0.005 up
f=0.005 down
f=0.01 up
f=0.01 down
a
0.6 0.7 0.8 0.9 1.0 1.1 1.2

0.0
0.2
0.4
0.6
0.8
1.0
1.2
Amplitude
f=0.001 up
f=0.001 down
f=0.005 up
f=0.005 down
f=0.01 up
f=0.01 down
b
0.6 0.7 0.8 0.9 1.0 1.1 1.2

0.0
0.2
0.4
0.6
0.8
1.0
1.2
Amplitude
f=0.001 up
f=0.001 down
f=0.005 up
f=0.005 down
f=0.01 up
f=0.01 down
c
101 103 105 107 109 1011
1/mi
0.000
0.005
0.010
0.015
0.020
0.025
????
f=0.001
f=0.005
f=0.01
zero line
d
Рис 3.1. Динамические характеристики системы: a — амплитудно-частотные характеристики 2-го
осциллятора при различных значениях f ; b — амплитудно-частотные характеристики 4-го осциллятора
при различных значениях f ; c — амплитудно-частотные характеристики 6-го осциллятора при
различных значениях f ; d — старший показатель Ляпунова Λ каждого осциллятора в цепочке при
различных значениях f
Fig 3.1. Dynamic properties of the system: a — frequency response of the 2nd oscillator for various values f ;
b — frequency response of the 4th oscillator for various values f ; c — frequency response of the 6th oscillator
for various values f ; d — the largest Lyapunov exponents Λ for every oscillator in the chain for various
values f
Как видно из сравнения рисунков, максимумы амплитуд соответствуют одинаковой частоте,
что означает синхронизацию колебаний осцилляторов. При низком значении f = 0.001 в
системе преобладает устойчивая динамика, однако с ростом амплитуды внешнего воздействия
появляются области неустойчивости. Например, с ростом значения f в окрестности частоты
Ω = 1 имеется хаотическая область, аналогичная той, что наблюдалась в [35; 36].
Выводы
В настоящей статье рассматривалась динамика связанных осцилляторов с убывающими в
геометрической прогрессии массами. Кроме того, предполагалось наличие не только нелинейных упругих, но и гистерезисных связей между осцилляторами. Исследовалась динамика указанной
системы в условиях периодического воздействия на первый элемент в цепочке.
Для анализа реакции системы на внешнее периодическое воздействие строилась амплитудно-
частотная характеристика для набора осцилляторов. Кроме того, численно строился энергетический спектр системы (зависимость сосредоточенной на n-м осцилляторе энергии от величины, обратно пропорциональной массе). Как следует из результатов вычислительных экспериментов, при достаточно низких значениях амплитуды имеет место ”почти” гармонический резонанс. При этом механическая энергия в основном локализуется на первом осцилляторе, а ее трансформация на остальные элементы цепи происходит по экспоненциальному закону. Однако при достаточно высоком значении амплитуды внешнего воздействия система реагирует хаотично.
В результате численных экспериментов были идентифицированы диапазоны частот внешнего воздействия, которым отвечает хаотическое поведение осцилляторов. При этом указанные отрезки существенно ´yже по сравнению с системой, в которой имеют место лишь упругие
взаимодействия. Отметим, что аналогичная ситуация наблюдалась для систем гистерезисно связанных осцилляторов Ван-дер-Поля [38]. Указанный эффект объясняется регуляризирующей
ролью гистерезисных звеньев, абсорбирующей часть энергии колебаний на каждом такте.
Из сравнения кривых энергетического спектра систем осцилляторов с гистерезисными
связями с системами осцилляторов с нелинейными упругими связями следует, что гипотеза Колмогорова находится в лучшем соответствии с данными, полученными в результате вычислительных экспериментов именно для систем с гистерезисными связями.

About the authors

M. E. Semenov

Voronezh State University; State Research Institute of Aircraft Systems (GosNIIAS)

Author for correspondence.
Email: mkl150@mail.ru
ORCID iD: 0000-0002-3361-1102

Doctor of Physical and Mathematical Sciences, professor of the Department of Digital Technologies

Russian Federation, 1, Universitetskaya Square, Voronezh, 394018, Russian Federation; 7, Viktorenko Street, Moscow, 125319, Russia Federation

A. V. Tolkachev

Voronezh State University; Voronezh State University of Forestry and Technologies named after G.F. Morozov

Email: tolkachev.akim@mail.ru
ORCID iD: 0000-0002-4968-5253

senior lecturer of the Department of Digital Technologies

Russian Federation, 1, Universitetskaya Square, Voronezh, 394018, Russian Federation; 8, Timiryazeva Street, Voronezh, 394613, Russian Federation

O. I. Kanishcheva

Voronezh State University

Email: oleka-olesya@mail.ru
ORCID iD: 0000-0002-1830-4091

Candidate of Physical and Mathematical Sciences, associate professor of the Department of Digital Technologies

Russian Federation, 1, Universitetskaya Square, Voronezh, 394018, Russian Federation

References

Supplementary files