О компьютерном моделировании обтекания высотных зданий и сооружений Текст научной статьи по специальности «Математика»

УДК 681.326

О КОМПЬЮТЕРНОМ МОДЕЛИРОВАНИИ ОБТЕКАНИЯ ВЫСОТНЫХ ЗДАНИЙ И СООРУЖЕНИЙ

ФИЛАТОВ Е.Ю., асп., ЯСИНСКИЙ Ф.Н., д-р физ.-мат. наук

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

Ключевые слова: математическая модель, ветровая нагрузка, методы вычислительной гидродинамики.

ON COMPUTER SIMULATION OF HIGH-RISE BUILDINGS AND CONSTRUCTIONS

FLOW OPERATION

E.Yu. FILATOV, postgraduate, F.N. YASSINSKIY, Ph.D.

The paper represents the mathematical model for calculating air motion relative to high-rise buildings and constructions. The authors calculate wind load over such constructions, using Sekundov's turbulence pattern, and try to find the solution on enclosed network variety.

Key words: mathematical model, wind load, calculating hydrodynamics approaches.

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

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

Рис. 1. Схема расположения узлов и здания

Гидродинамические уравнения и разностные схемы.

Решение задачи состоит из следующих этапов:

1. Вычисление поля скорости. Для вычисления поля скорости используется система ю-у (вихрь -функция тока). Расчет ведется следующим образом. Имея начальное распределение поля вихря ю, с помощью уравнения Пуассона находим поле функции тока у:

2 2 д у д у —2 + —2 = -ю. дх2 ду2

Далее находится поле скоростей и (по определению функции тока):

их ; П =-дУ.

х ду у дх

По полю скоростей ищется поле вихря ю через уравнение Гельмгольца:

дю

-U,

дю дх

,, дю Uy= v

у ду

д ю

2

д ю

чдх 2 ду \

Рассмотрим теперь более подробно описанные шаги и представим разностные схемы к приведенным уравнениям.

В данной реализации используется многосеточный метод, который будет описан ниже. За основную берется сетка размером 62х28 узлов (плюс узлы четырех границ). Внутрь этой области помещается здание размером 5х5 (больших) шагов по расстоянию. Через Л будем обозначать шаг по расстоянию, через т - шаг по времени.

Начальные условия. Скорость и и вихрь ю принимаются равными 0 по всей области (за исключением границ). Функция тока у изменяется линейно по вертикали от наименьшего до наибольшего значения (см. рис. 1).

Гоаничные условия. Для у на правой границе используется граничное условие 2 рода (первая производная равна 0), на остальных границах -1 рода (со значениями, заданными в начальных условиях). На границе здания у становится постоянной, равной полусумме максимального и минимального значений у. Для ю на границах используется граничное условие 2 рода, на границе со зданием -1 рода. Скорость на границах принимается равной заданной в условии скорости воздушного потока, а на границе со зданием - равной 0.

Уравнение Пуассона для нахождения функции тока у решается с использованием метода верхней релаксации:

уновое = (1 - 0)у,-, у + 4 (у/+1,1 + +у/-1,1 + у/, 1+1 + У/, 1 -1 + h2ю),

где 0 - некоторая константа.

Поле у обновляется до тех пор, пока не станет с некоторой погрешностью стационарным.

Для нахождения поля скорости и используются центральные производные по расстоянию:

у/,1+1 - у/, 1 -1 ,, = у/-1,1 - у/+1,1 , иу/, 1 =

Uxi, j =-

2h

2h

Уравнение Гельмгольца для нахождения поля вихря ю решается в неявной схеме с использованием метода расщепления. Шаг по времени делится на две части. В первой части расчет ведется в направлении Ох, во второй - в направлении Оу.

В первой части шага решается уравнение

дю ,, дю д2ю — + Ux — = V—г-. дt x дx дx2

Разностная схема для этого уравнения выглядит так:

k+1 \ j

и/, 1

( д2 +1 д ю 2

= -и

xi, 1

к+1

| д2 +2 д ю 2

дх2

дх 2

к+1 2

/+1,1

- 2ю

к+1 /. 1

к+1 э/-12

Ь2

Для улучшения устойчивости используем про-тивоточные производные:

Цх',< (I) =

и

-1,1

х/, 1

Ь

их

/+1,1

/. 1

Ух-,, 1 > 0,

Ух1,1 < о.

х/,' Ь

Во второй части шага решается уравнение

., дю д2ю Уу — = V—2. у ду ду2

дю

"дГ

Разностная схема для этого уравнения строится аналогично предыдущему.

Поскольку используется неявная схема, то для разрешения полученных разностных уравнений используется метод прогонки. Разностные уравнения для первой части шага приводятся к виду к+1 , к+2 к+2 аю/-12>-+ь,ю,, 12 +сю/+12>- = Ъ■

Полагая

к+1

к+

ю/-12у = -¡ю, у 2 + М/, находятся коэффициенты — и М :

— = с/ ■ М = ъ- аМ /

-/+1 =--;-г" ■ М/+1 =

а- +ь.

а- +ь.

к+1 1

-¡, М, и 1,1 на границах находятся исходя из граничных условий:

1) для граничных условий 1 рода -

— = 0;

М1 = ю0,1 ■

к+1 _ к юЫх+1,1 = юЫх +1,1,

где Ых - число узлов в направлении Ох; 2) для граничных условий 2 рода -

- = 1; М1 = 0;

к+2 °Ых +1,1

М

+1

Ы +1

-1

В результате находится новое значение поля вихря ю.

2. Вычисление турбулентной вязкости по Се-кундову. Рассматриваемое движение турбулентно, поэтому необходимо рассчитывать также поле турбулентной вязкости. В данной реализации для расчета турбулентной вязкости использовалась модель Се-кундова. В этой модели поле турбулентной вязкости п описывается следующим уравнением:

дv ,, ду ,, ду д | /

— + Ух — + иу — = —I (у

дt х дх у ду дх

ч ду А — ' дх

д | ( +— I (

ду Г

ч ду —

' ду

¿-I о-Я-(у

—гтш

-ру)у,

где у/ - ламинарная вязкость; 1_тт - кратчайшее расстояние до ближайшей неподвижной твердой поверх-

ности;

О =

Ъ (г) = 0,2

а = 2; р = 0,06; у = 50.

I2 -1,47г +1

При составлении разностной схемы для второго и третьего слагаемых левой части уравнения используются противоточные производные. Разностная схема для первых двух слагаемых правой части выглядит следующим образом:

у / (у+1- 2у/+у 1-1)+а (у'+12 - 2у ¡2+у,-12)).

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

На левой границе вязкость принимается равной ламинарной, на остальных используется граничное условие 2 рода (первая производная равна 0). На границе со зданием вязкость принимается равной 0.

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

3. Использование сгущающейся сетки и асинхронного интегрирования. В разных местах вокруг здания изменения полей происходят неравномерно. Поэтому было предложено уменьшить шаг по времени и шаг по расстоянию там, где изменения происходят более динамично. Для этого вся рассматриваемая область была разделена на 3 зоны, как показано на рис. 1. Первая -самая широкая - охватывает вторую, которая, в свою очередь, охватывает третью, содержащую здание. Такое расположение зон было выбрано исходя из того, что наиболее быстрые и важные изменения происходят ближе к зданию. Если шаг по расстоянию первой зоны равен Ь, а по времени - т, то во второй зоне шаги были выбраны соответственно Ь/2 и т/2, а в третьей - Ь/4 и т/4. Вариацию шага по расстоянию мы называем методом сгущающейся сетки, а по времени - асинхронностью интегрирования. Назовем шаг Ь большим шагом по расстоянию, а шаг т - большим шагом по времени.

Распределение узлов на границе сеток представлено на рис. 2.

N»-1

N,6

М1Ь+1

0 12 3 4 5

N,3-1 N,3+1

Рис. 2. Схема распределения узлов для вложенных сеток

Рис. 3. Полученные поля скорости, линий тока и турбулентной вязкости

Расчет в первой зоне производится в узлах до Л1а-го и до Л/1Ь-го. Значения в узлах (Л/1а+1)-го столбца и (Л1Ь+1)-й строки принимаются от второй зоны. Во второй зоне расчет ведется начиная с 1-й строки и 1-го столбца. Значения в 0-й строке и 0-м столбце принимаются с предыдущей зоны, используя линейную интерполяцию.

Асинхронность интегрирования выполнена следующим способом. В первой зоне делается 1 (большой) шаг по времени т, во второй - 2 шага по т/2, в третьей - 4 шага по т/4, а затем происходит обмен значениями между зонами. Иными словами, обмен значениями между зонами происходит через каждый большой шаг по времени.

4. Расчет нагрузки. Чтобы найти ветровую нагрузку, действующую на здание, необходимо найти распределение давления вокруг здания. Для этого решается уравнение Пуассона

д2Р

д2Р

(

дх2 ду2

■ = -р

5Цх

дх

дЦу +№

ду дх I ду

2 ^

где р - плотность воздуха.

Уравнение решается подобно описанному выше уравнению Пуассона для нахождения функции тока у.

По полученному полю давления вокруг здания находится удельная нагрузка на здание (т.е. сила воздушного потока на единицу высоты здания) как разность давлений на левой Г/ и правой Гг границах здания:

д2

| Рду - | Pdу.

Г/ Гг

Интегралы вычисляются методом трапеций.

Результат работы программы представлен на рис. 3. Кривыми изображены линии тока, стрелками -вектора скорости, темным выделена область высокой турбулентной вязкости. Воздушный поток огибает здание, создавая за зданием область турбулентных пульсаций. Нагрузка на здание была рассчитана и получилась равной 7573 Н/м.

Список литературы

1. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980.

2. Численные методы и параллельные вычисления для задач механики жидкости, газа и плазмы / Э.Ф. Балаев, Н.В. Нуждин, В.В. Пекунов и др.: Учеб. пособие. -Иваново, 2003.

3. Лойцянский Л.Г. Механика жидкости и газа. - М.: Наука, 1970.

Филатов Евгений Юрьевич,

ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», аспирант кафедры высокопроизводительных вычислительных систем, телефон (4932) 26-98-29.

Ясинский Федор Николаевич,

ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,

доктор физико-математических наук, профессор, зав. кафедрой высокопроизводительных вычислительных систем,

телефон (4932) 26-98-29.