Моделирование пространственного поля скоростей атмосферной турбулентности приземного слоя атмосферы Текст научной статьи по специальности «Физика»

__________УЧЕНЫЕ ЗАПИСКИ ЦАГИ

Том XXX 199 9

М3—4

УДК 551.551

МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННОГО ПОЛЯ СКОРОСТЕЙ АТМОСФЕРНОЙ ТУРБУЛЕНТНОСТИ ПРИЗЕМНОГО СЛОЯ АТМОСФЕРЫ

В. Е. Грязин

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

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

Методы, используемые для моделирования атмосферной турбулентности, зависят от особенностей рассматриваемой задачи. Они определяются режимом полета летательного аппарата, его конструкцией и геометрическими размерами. Если геометрические размеры летательного аппарата малы по сравнению с характерным размером турбулентных вихрей, то можно считать, что во всех точках его поверхности скорости порывов ветра одинаковые. Если при этом траектория полета летательного аппарата близка к прямолинейной, задача моделирования атмосферной турбулентности может быть сведена к моделированию одномерных случайных процессов, для описания которых используются их спектральные плотности. Такие методы моделирования просты, требуют малых затрат машинного времени.

В том случае, когда геометрические размеры летательного аппарата соизмеримы с масштабом турбулентности и/или траектория полета летательного аппарата существенно отличается от прямолинейного полета, появляется необходимость моделирования атмосферной турбулентности

в виде трехмерного случайного поля. Это позволяет учесть дополнительные приращения аэродинамических сил и моментов, возникающих за счет неравномерного распределения порывов ветра по поверхности летательного аппарата.

1. Основные предположения о характере атмосферной турбулентности, используемые при разработке алгоритмов моделирования.

В общем случае трехмерный вектор турбулентности \fVxg, №уё, является функцией трех пространственных координат Х£,у8,г8 и времени ?.

Обычно при анализе результатов метеонаблюдений используется земная система координат, ориентированная по направлению среднего ветра. Таким образом, ось направлена по направлению среднего ветра, а уё и

2^ являются соответственно вертикальной и поперечной осями. Направление среднего ветра в нижнем слое высот считается неизменным, что характерно для нейтрального состояния атмосферы, встречающегося приблизительно в 70% метеонаблюдений.

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

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

Случайное поле называется изотропным, если корреляционные функции зависят только от расстояния между двумя точками и не зависят от положения пары точек. Предположение об изотропности поля атмосферной турбулентности справедливо только на больших высотах, где не сказывается влияние земной поверхности. По данным, приведенным в работе [2], поле атмосферной турбулентности становится изотропным на выЛ 10

сотах, превышающих Н = ЮОО/г0’ м, где /?о М — коэффициент шероховатости подстилающей поверхности.

Гипотеза о «замороженном» поле атмосферной турбулентности предполагает, что за время движения рассматриваемого объекта сквозь поле атмосферной турбулентности оно не изменяется по времени. Эта гипотеза справедлива при относительно быстром движении объекта, когда турбулентные вихри не успевают разрушаться. В работе [1] приведена оценка времени распада турбулентных вихрей тр «100 с. Поэтому, если объект

проходит характерный размер вихрей за интервал времени значительно меньший, чем Хр, вихри можно считать «замороженными». В работе [3]

приводится следующее условие применимости этой гипотезы:

у/и >1 ь/300, где ‘V —скорость движения летательного аппарата относительно воздушного потока в /-м направлении; I/ — скорость среднего

ветра; 1Ь — характерный размер вихрей атмосферной турбулентности (масштаб) в г'-м направлении в метрах. Учитывая вышеизложенное, можно предположить, что для большинства случаев, представляющих практический интерес для авиации, гипотеза о «замороженном» поле атмосферной турбулентности применима. Использование этой гипотезы позволяет исключить из рассмотрения время и уменьшить число аргументов поля атмосферной турбулентности до трех пространственных координат.

2. Алгоритмы моделирования гауссовского изотропного поля атмосферной турбулентности. Представление атмосферной турбулентности в виде гауссовского изотропного поля значительно упрощает задачу моделирования. В этом случае поле атмосферной турбулентности полностью описывается двумя характеристиками: математическим ожиданием и матричной корреляционной функцией Я [р;/]. Причем элементы матрицы Я

г — вектор, соединяющий две точки пространства | г

гп rj — проекции вектора г на оси , у;

8,у — символ Кронекера, 5у = 1 при / = у, 5^ = 0 при / * у ;

■аг-, а у — среднеквадратичные значения компонентов атмосферной турбулентности, в дальнейшем для упрощения принимается аг , а у = 1.

определяются двумя пространственными корреляционными функциями: продольной/(г) и поперечной £(г) [4]:

і, у -1, 2,3,

где

Вид функций /(г) и g(r) для модели атмосферной турбулентности в форме Кармана приведен в работе [3]. Используя преобразование Фурье, по приведенной матричной корреляционной функции можно получить следующее выражение для элементов матрицы спектральной плотности в форме Кармана [4]:

55 Г025;, - .•'] .

1---------^---2-Г7Т6-’ (2)

1 36 п2а(аЬ)2В [(о/,)-2 +02]17/6

где

О — пространственная частота, О = (Пь С22, П3)1;

Ь — интегральный масштаб турбулентности; а - 1,339 — постоянная Кармана.

Для гауссовского изотропного трехмерного поля атмосферной турбулентности можно синтезировать векторную параметрическую модель в виде [5]:

N . ,.

Щх) = (т8(х), Щ(х), ^(х))т= °(Д,))’ (3)

^ к=1

где

£(*, С1^) = П^)> £,2(Х’СЪ^Ъ> --- трехмерное

случайное поле, зависящее от аргумента х = [хё, у8,г8)'т и случайного

трехмерного вектора ; N— число независимых реализаций случайного поля.

Если векторы статистически независимы и одинаково распределены с плотностью распределения (р(О), а поле %(х, обладает

свойствами

М

о, м

Ру-’ (4)

то математическое ожидание и матричная корреляционная функция векторного поля атмосферной турбулентности и его параметрической модели (3) совпадают при всех Ы, а конечномерные распределения IV(х) при N -» оо нормальны. Поэтому модель (3) при достаточно больших N может быть использована для моделирования гауссовского векторного поля атмосферной турбулентности. В качестве используются про-

странственные гармоники со случайными частотой О и фазой ф, равномерно распределенной в интервале [0,2л]. Амплитуда гармоник (3 зависит от частоты О и случайного параметра г.

4/(х, 0) = рДг, □)8Іп[ 0Т*-і-ф], / = 1,2,3,

где

т

Од: — скалярное произведение;

г = (2], ^2, 23) — независимые случайные величины, м(гу )= О,

О = (□], 02, ^з)Т — трехмерный случайный вектор, независимый от г с плотностью распределения ф(О) > 0 .

Для выполнения условий (4) элементы матрицы А{а^) определяются из выражения

5(0) = Л(0)Ф(0М(0)Т.

(6)

Здесь 5(0) — матричная спектральная плотность поля атмосферной турбулентности, Ф(О) — диагональная матрица, все элементы которой равны произвольной функции ф(О), обладающей свойствами:

По аналогии с моделированием случайных процессов матрица А называется частотной характеристикой пространственного формирующего фильтра.

Модель случайного векторного поля (3), (5) содержит две неизвестные функции: А(0.) и ф(О). Для упрощения моделирующих алгоритмов произвольная функция ф(О) принимается равной спектральной плотности одной из компонент моделируемого поля, например: ф(О) = 5] 1(0).

Из свойств спектральной плотности следует, что такая функция при

о, = 1 удовлетворяет требованиям, предъявляемым к функции ф(О). После выбора функции ф(О) можно найти соответствующие ей элементы матриц А{О), используя выражение (6). В работе [5] показано, что при предположении о корреляции всех трех компонент турбулентности друг с другом матрица А{€1) является треугольной с элементами:

аи - и аі\ = -

О ] О 2 пі+пі’

а 22

іпцп31

О 2 ^ 3

О і о

1**3

+ П§

-8І8п(П3)|£2|П2 . _ л

а32 ------------------, а33 - 0.

02+^3

Таким образом, для моделирования изотропного случайного поля атмосферной турбулентности по формулам (3), (5) достаточно иметь алгоритмы, моделирующие случайные величины 2 и случайный вектор О . Как отмечалось выше, независимые величины г, должны удовлетворять условиям м\zj )= 0, = 1. Этим условиям, например, удовлетворяет слу-

чайная величина, равновероятно принимающая значения ±1.

Плотность распределения вектора ф(О) = 5] ] (О) имеет вид:

ср(С2і, 02, Оз) =

551

1

36 %2а(аЬ)2/3

{аьу1 + а2

17/6 '

(8)

Обозначим Х = ^02 +О3 — модуль двумерного вектора (о2, П3), распределенного, как видно из выражения (8), изотропно. Интегрируя ф(£2) по О] е [-оо, со], находим совместную плотность распределения 02 и О3 :

АХ1

9ті(а£)2/3і\2 + (аІ)_2У/3

(9)

Вектор (о2>^з) с плотностью распределения (9) в полярных координатах моделируется по формуле:

(£22, Оз)Т=Яе.

Здесь е = (е1,е2)Т — реализация изотропного направления на плоскости, которое задается полярным углом ці, равномерно распределенным на интервале [0,2л]. Моделирующий алгоритм имеет вид:

е\ - соз(\|/); е2 = 5т(у).

Плотность распределения модуля л. —> /о(Х) в полярных координатах имеет вид:

М1) = 8кт = 2пХДЪ) =

8

9(^)2/3(я.2+(111 Г2)7'3’

(10)

где 5^ — длина окружности радиуса 1. Для получения алгоритма моделирования случайной величины с плотностью распределения (10) используется метод обращения [5]. Для этого функция распределения величины X приравнивается равномерно распределенной величине у!:

/?о(0 = //о(^)^=^{^)2+1]-4[[^)2+1]”1/3+з} = у1. (П)

о

( 2 V1/3

Производя замену переменной у = (гаЬ) +1 , получаем неполное

уравнение четвертой степени для нахождения величины у :

F(у) = у* -4у + 3(l - У] )= 0 (12)

Искомое решение этого уравнения при изменении у] в интервале [0, 1] также лежит в интервале [0, 1]. Решение уравнения (12) можно найти численно, используя итерационный метод последовательных приближений.

Далее по формуле г = (1/аЬ)^\/у3 - lj находится величина модуля X, и

величины Q2 и определятся по формулам:

Q2 = г cos (2тгу2), fi3 — /- sin (2тху2), ' (13)

где у 2 =Rav[0,1] — равномерно распределенная величина и независимая от у].

Первая компонента вектора Q -» Qj моделируется как случайная величина с плотностью распределения:

z X 9^02^ _55______________________1______

/,(l) /w

л 2 —2 ‘

Здесь b -X +{aL) . Анализ выражения (14) показывает, что случайная

величина Qj может быть представлена в виде Qj = br\, где г| — случайная

величина с плотностью распределения:

(15)

(i+л J

Для нахождения моделирующего алгоритма величины г| используется метод обращения. Функция распределения величины г| приравнивается равномерно распределенной величине у3:

После подстановки ті = г§(ц); й?г| = ф/соэ2 ц; ц є [-я/2, я/2] получаем уравнение:

(16)

Интеграл в уравнении (16) в элементарных функциях не выражается, поэтому значение верхнего предела [3, соответствующее выпавшему значению случайной величины у3, определяется численно. После этого случайная величина О] вычисляется по формуле:

В результате можно сформулировать следующую последовательность моделирования изотропного поля атмосферной турбулентности.

а) при выпавшем значении величины У1 = Яау [0,1] итерационным методом решается алгебраическое уравнение четвертой степени (12) и вы-

б) по формулам (13) вычисляются значения второй и третьей компоненты вектора О;

в) «выбрасывается» значение случайной величины у3 и. по формулам (16), (17) вычисляется значение первой компоненты вектора Г2.

2. Вычисляются случайная амплитуда и фаза гармоники:

а) по формулам (7) вычисляются коэффициенты треугольной матрицы А\

б) «выбрасываются» значения случайных величин 2} и по формуле

(5) вычисляется случайная амплитуда гармоники;

в) «выбрасывается» случайная фаза гармоники ср = Лау [0,2л].

3. Пункты 1 и 2 повторяются N раз, и вычисленные значения запоминаются. :

4. При фиксированном значении пространственного положения точки х = (х^, угпо формулам (3), (5) вычисляются компоненты случайного вектора атмосферной турбулентности.

5. Пункт 4 повторяется многократно для различных значений пространственного положения рассматриваемой точки или для нескольких точек летательного аппарата.

(17)

1. Вычисляется случайный вектор О = (Г2|, С12, Оз )т :

числяется модуль двумерного изотропного вектора (о2, Пз);

3. Учет анизотропии поля атмосферной турбулентности в нижнем слое высот. В разделе 2 описаны алгоритмы моделирования изотропного поля атмосферной турбулентности. Эта модель пригодна только для больших высот. В работе [2] отмечалось, что поле атмосферной турбулентно-

которрй определяется интегральный масштаб;

р„ (А/, 0) — автокорреляционная функция г-й компоненты атмосферной турбулентности при нулевом временном запаздывании и расстоянии между двумя точками измерения А/ = Ах^, Ау^, Аг^.

Учет анизотропии поля атмосферной турбулентности осуществляется путем введения зависимости масштабов каждой компоненты атмосферной

турбулентности ] Ц от высоты полета [3].

Моделирование анизотропного поля атмосферной турбулентности, имеющего в приземном слое атмосферы эллиптическую анизотропию, осуществляется путем сжатия, полученного в разделе 2, изотропного поля по осям пространственных координат х8,уё,гё. Для этого в выражении (10) в качестве исходного берется продольный интегральный масштаб турбулентности, равный 280 м. При Ь = 280 м вычисляются случайные амплитуды, частоты и фазы гармоник по пунктам 1—3. Далее по формулам (3), (5) вычисляются значения компонент скорости ветра в исходной точке пространства 2^. Обозначим их соответственно

. При перемещении рассматриваемой точки из исходного

положения в новую точку пространства х^, у^ в анизотропном поле атмосферной турбулентности происходит изменение положения точки на величину Н = -г^). В изотропном поле, масштабы ко-

торого отличаются от масштабов анизотропного поля, эквивалентное перемещение точки будет ]1, =Н 28о/-7^ . Другими словами, для того чтобы получить в изотропном поле с большим, чем в анизотропном поле, масштабом /-й компоненты турбулентности такое же (в статистическом смысле) изменение скорости ветра, необходимо, соответственно, пройти

А 1 О

сти изотропно только на высотах Н >1000/^’ м. Масштабы компонент атмосферной турбулентности равны:

280 м при / = у; 140 м при /* /;

Здесь

— направление, вдоль которого определяются масштабы турбулентности;

— компонента атмосферной турбулентности, для

в у'-м направлении большее расстояние. Поскольку масштабы являются непрерывными функциями высоты, то необходимо сжимать изотропное поле на каждой высоте полета со своим коэффициентом сжатия:

где Vj(t) — скорость перемещения точки в у-м направлении; / — время

движения рассматриваемой точки. Необходимо также учитывать, что за время t происходит не только перемещение рассматриваемой точки в системе координат уё, , фиксированной на поверхности земли, но и перемещение самого поля атмосферной турбулентности, которое переносится средним ветром и. Поэтому суммарное перемещение рассматриваемой точки относительно поля атмосферной турбулентности в продольном направлении составляет:

С учетом вышеизложенного окончательное выражение для моделирования г-й компоненты анизотропного поля атмосферной турбулентности принимает вид:

Здесь а,- — среднеквадратичные значения /-й компоненты атмосферной турбулентности, зависящие от высоты и скорости среднего ветра [4], Ве-

пример нулевые, значения. На каждом шаге интегрирования уравнений

затем вычислять компоненты атмосферной турбулентности .

4. Результаты математического моделирования. Для проверки совпадения статистических характеристик моделируемого поля атмосферной турбулентности с результатами метеонаблюдений были проведены статистические расчеты коэффициентов корреляции компонент атмосферной турбулентности. На рис. 1 приведены коэффициенты корреляции про-

280/4 . В случае движения точки сквозь поле атмосферной турбулентности эквивалентное перемещение необходимо вычислять путем интегрирования:

личины в начальный момент времени могут иметь произвольные, на-

движения будет необходимо вычислять масштабы ]Ь( перемещения Ні и

Н, М

200

150

100

50

у 0,2 0,4 0,6 0,8 р

Рис. 1

дольной и поперечной составляющих атмосферной турбулентности в двух, разнесенных по высоте на 25 м, точках в функции высоты расположения верхней точки. Приведенные результаты получены путем осреднения по 20 000 реализациям поля атмосферной турбулентности. Для сравнения на рисунке приведены аналогичные коэффициенты корреляции, полученные по результатам метеонаблюдений на метеовышке Института экспериментальной метеорологии (ИЭМ).

На рис. 2 приведены коэффициенты корреляции про- Л дольной и поперечной составляющих атмосферной турбу- 08 лентности в функции приведенного продольного рас- 0,б

стояния между двумя точками

г = 0,747 Ц ) (тем- 0,4

ные метки). Результаты полу- 0,2

чены путем осреднения по 20 000 реализациям поля ат- о

мосферной турбулентности.

Для сравнения на рисунке I 0 2 приведены исходные зависимости продольной /(г) и поперечной g(r) корреляционных функций [3]. На рис. 2 также приведены коэффициенты корреляции продольной и вертикальной компонент атмосферной турбулентности в функции г для двух точек, расположенных на консолях крыла самолета с размахом / = 30 м (светлые метки). Эти коэффициенты рассчитаны при моделировании длительного прямолинейного полета самолета сквозь одну реализацию поля атмосферной турбулентности на двух высотах (Н= 100; 50 м).

Как видно из приведенных результатов, коэффициенты корреляции моделируемого случайного поля атмосферной турбулентности хорошо

)» £(/■) • *. Расчет по модели

о о

Рис. 2

совпадают с исходными. Однако совпадение коэффициентов корреляции, полученных по одной реализации, удается получить только при большом числе гармоник, используемых для моделирования атмосферной турбулентности N>50.

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

Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 98-01-000174).

ЛИТЕРАТУРА

1. PanofskyY. A.,DittonJ. A. Atmospheric turbulence. Models and methods for engineering applications.— John Wiley & Sons, New York and oth.—

1984. Vol. XIX.

2.: Characteristics of atmospheric turbulence near the ground//Engineering ,

Sciences Data Item Number 74030. ESDU-1974. .

3. Characteristics of atmospheric turbulence near the ground//Engineering Sciences Data Item Number 75001. ESDU-1975.

4. Etkin B. Turbulent wind and its effect on flight//!. Aircraft.— 1981.

Vol. 18, N5. . r .

5. Шалыгин А. С., Палагин Ю. И. Прикладные методы статистического моделирования.— М.: Машиностроение.— 1984.

Рукопись поступила 3/IX 1998 г.