CC BY
10
УДК 681.326
ПРЕИМУЩЕСТВО ИСПОЛЬЗОВАНИЯ МЕТОДОВ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ ПРИ РАСЧЕТЕ ВЕТРОВОЙ НАГРУЗКИ, ДЕЙСТВУЮЩЕЙ НА СООРУЖЕНИЯ РАЗЛИЧНОЙ ФОРМЫ И НА СТРОИТЕЛЬНЫЕ КОМПЛЕКСЫ
ФИЛАТОВ Е.Ю., асп.
Решается задача расчета ветровой нагрузки, действующей на здания и сооружения, методами вычислительной гидродинамики. Предлагаются методы ускорения вычислений и повышения их точности. Показываются недостатки использования СНиП при решении задач расчета ветровой нагрузки и преимущество изложенных методов.
Ключевые слова: численное моделирование, поле скорости, граничные условия, модель турбулентности, ламинарный подслой, метод сгущающихся сеток.
ADVANTAGES OF USING THE METHODS OF COMPUTATIONAL FLUID DYNAMICS WHILE ACCOUNTING THE WIND LOAD ON CONSTRUCTIONS OF DIFFERENT FORMS AND BUILDINGS COMPLEXES
E.Yu. FILATOV, Post Graduate Student
The author investigates the accounting problem of wind load on buildings by means of fluid dynamics methods. The article represents computational speedup methods and accuracy increase methods. The disadvantages of SNiP usage are represented while solving wild load accounting problems. The article shows the advantages of the introduced methods.
Key words: numerical simulation, speed field, boundary conditions, turbulence model, laminar sublayer, method of infill nets.
Введение. Для расчета ветровой нагрузки, действующей на здания и сооружения, традиционно применяют аэродинамические формулы, изложенные в Строительных Нормах и Правилах (СНиП) [3]. В случаях, когда решение задачи с помощью формул СНиП затруднено, проводят натурный эксперимент - продув уменьшенной модели сооружения или строительного комплекса в аэродинамической трубе. Такой эксперимент дает достаточно точные результаты, однако он дорогостоящ и имеет ограничения как на размеры продуваемых объектов, так и на скорость потока, а следовательно, и на число Рейнольдса. Что касается использования формул СНиП, то они дают результаты невысокой точности, в таблице СНиП присутствует ограниченное множество форм сооружений, а также в формулах СНиП очень грубо учитывается окружающая застройка. Ниже будут приведены примеры, в которых использование формул СНиП дает заведомо неверный результат. Все это говорит в пользу иного метода расчета ветровой нагрузки, обширное применение которого стало возможным благодаря экспоненциальному росту мощностей вычислительной техники. Это численное моделирование аэродинамических процессов на основе методов вычислительной гидродинамики.
Решение задачи. Далее будет решаться следующая задача. Дано одиночное здание в форме параллелепипеда и дан профиль скорости ветра, направленный параллельно одной из сторон здания. Требуется рассчитать поля скорости и давления вокруг данного здания и определить силу давления воздуха на здание. Обобщение данной задачи на несколько зданий и на здания и сооружения различной формы не представляет особого труда.
Данная задача решается в системе «давление -скорость» [1, 2]. Для нахождения поля скорости
и решаются уравнения Навье-Стокса:
j=1 dYj
v>
d2U,
By,2 p dx
(1)
у=1 ^7
Для нахождения поля давления Р решается уравнение Пуассона:
AP = _р
— + > U,---------------------------v>
Pit ^ ■* Яу . ^
dD
d2D
3 3
dt 1=1
DU, dUj;
dYj
j=1 dYj
Q
(2)
(3)
0 = > >
м у=1 дху дх1
0 = ^ и = > 1т (4)
I=1 дХ1
Здесь р - плотность; V = Vmol + vturb - кинематическая вязкость (сумма молекулярной и турбулентной вязкости); О - поле дивергенции скорости.
В качестве начальных условий для поля давления на всей рассчитываемой области берется некоторое начальное значение Р0, а для поля скорости - значение, соответствующее профилю скорости. Профиль скорости вычисляется по формуле
H,
ref
(5)
где vz - скорость ветра на высоте z; vref - скорость ветра на высоте рассматриваемой области Href; а = 0,26 - показатель степени.
В качестве граничных условий для поля скорости на выходной границе используется граничное условие II рода, при котором производная по нормали к границе равна 0. На остальных свободных границах используются граничные условия
I рода, при которых скорость соответствует заданному профилю, а на здании используется граничное условие I рода, при котором скорость равна 0. Для поля давления на выходной границе ставится мягкое граничное условие, при котором вторая производная по нормали к границе равна 0. На остальных границах (на свободных границах и на здании) используются граничные условия II рода, при которых первая производная по нормали к границе равна 0.
Поскольку рассматриваемое течение турбулентно, требуется рассчитывать также поле турбулентной вязкости. Для этого используется модель турбулентности А.Н. Секундова [2, п.1.2.2]. Она позволяет получить приемлемые результаты при невысоких вычислительных затратах. Турбулентная вязкость vturb в данной модели рассчитывается из уравнения
ду
дt
1=1 дх1
=>— (
г)у . ^
1=1 дх1
vmol
ЧигЬ)
игЬ
дХ:
ЧигЬ
8у
то!
turb’
О =
1,1
ди1
дх.
ди
_____}_
дхI
1/2
(6)
(7)
/У
г2 +1,47 г + 0,2.
1
(8)
(9)
г2 -1,47га = 2, р = 0,06, у = 1.
Здесь О - величина, оценивающая интенсивность деформации среды; 1-тП - кратчайшее расстояние до ближайшей твердой поверхности; 1 (-)- эмпирическая функция; а, р, у - эмпирические константы.
Для того чтобы поле турбулентной вязкости сформировалось, требуется задать его начальное значение отличным от 0, например
уЮгЬ° = 1000у то!. (10)
В качестве граничных условий на входной границе, на земле и на здании используются граничные условия I рода, при которых турбулентная вязкость равна 0, на остальных свободных границах используются граничные условия II рода, при которых производная по нормали к границе равна 0.
Отдельно проводится расчет ламинарного подслоя пограничного слоя на поверхности здания, параллельной направлению ветра [2, п. 1.3]. Толщина ламинарного подслоя значительно меньше шага сетки, и ею можно пренебречь. Поэтому считается, что граничные узлы сетки попадают не на боковую поверхность здания, а на границу ламинарного подслоя, расположенного на боковой поверхности здания. При этом скорость на граничных узлах сетки равна скорости на границе ламинарного подслоя. Таким образом, при расчете полей сначала аналитически находится скорость на границе ламинарного подслоя, а затем она используется в качестве граничного условия I рода при расчете поля скорости.
Сила, действующая на здание, находится как сумма силы давления воздуха на здание и силы трения воздуха о здание:
Г = Гр
(11)
Сила давления вычисляется как интеграл давления по поверхности здания:
Гр = -.. Pnds, (12)
г
где Б - поверхность здания; П - вектор нормали к элементарной поверхности ds, направленный от здания.
Сила трения вычисляется по формуле
ди
дп
-ds.
(13)
Расчеты показывают, что сила трения много меньше силы давления (в одной из задач сила трения составляет 3,9-10-4 от силы давления), поэтому ею можно пренебречь.
Методы оптимизации вычислительного процесса. Изменение искомых полей происходит неравномерно во всей рассматриваемой области. Наиболее динамично изменение происходит вблизи здания. В удалении же от здания поля меняются довольно медленно. Кроме того, вблизи здания желательно иметь большую точность расчета, т.к. поля вблизи здания в значительной степени влияют на расчет полей во всей области. Исходя из этого, расчет полей был бы более эффективным, если бы шаг сетки и шаг по времени вблизи здания были меньше, чем в удалении от здания. Этого позволяют добиться методы сгущающихся сеток и асинхронного интегрирования по времени.
Метод сгущающихся сеток состоит в том, чтобы в разных частях рассматриваемой области выбирать разный шаг по расстоянию. Так, вся область делится на 3 зоны (рис. 1). Малая зона расположена вокруг здания, средняя зона - вокруг малой и большая - вокруг средней. В большой зоне шаг сетки равен И, в средней - И/2, в малой - И/4.
Рис. 1. Метод сгущающихся сеток и метод асинхронного интегрирования по времени
Кроме того, при уменьшении шага по расстоянию для сохранения устойчивости требуется также уменьшать шаг по времени, причем примерно в той же степени, что и шаг по расстоянию. Так, если в большой зоне шаг по времени равен т, то в средней зоне шаг по времени равен т/2, в малой зоне - т/4. Метод, при котором в разных частях расчетной области используются разные шаги по времени, называется методом асинхронного интегрирования по времени. Шаг по времени т называется большим шагом по времени.
Упрощенно алгоритм расчета полей в условиях сгущающихся сеток имеет следующий вид:
1. Расчет большого шага по времени в каждой зоне.
2. Обмен значениями полей на границах между зонами.
Под расчетом большого шага по времени понимается расчет 4-х шагов по времени т/4 в малой зоне и расчет 2-х шагов по времени т/2 в средней зоне. При расчете полей на границах зон используется граничное условие I рода и значения, полученные от соседней зоны.
Существуют два подхода в методе сгущающихся сеток. Они различаются распределением узлов на границе между зонами. При 1-м подходе (рис. 2) граничные (0-е) узлы меньшей (д + 1)-й зоны не совпадают с узлами большей д-й зоны. При 2-м подходе (рис. 3) граничные (0-е) узлы меньшей (д + 1)-й зоны совпадают с (Л - 1)-ми узлами большей д-й зоны. Здесь граничными узлами большей д-й зоны являются А-е узлы.
А+2-А+1 ■ А А-1 •
5
4 ■"і
О 2
0
д+1
Рис. 2. Метод сгущающихся сеток: 1-й подход
А+1
А
А-1
А-2
4
3 О
1
0
Рис. 3. Метод сгущающихся сеток: 2-й подход
Таким образом, обмен значениями полей на границах схематично можно изобразить следующим образом.
Для 1-го подхода:
[д ]А = в[д+1]1,
3[Я+1]0
Лд]А-1
Ля ]А
2
(14)
(15)
Ля ]А
= в[
[Я+1]2,
®[д+1]0 = ®[д ]А-1.
(16)
(17)
Для 2-го подхода: д +1 ^ д : <
д ^ д +1:
Здесь Э - данное поле; в квадратных скобках указывается номер зоны; последний индекс - номер узла по одной из координат.
Если сравнивать эти два подхода, то реализация 2-го подхода будет заметно проще. Это связано с тем, что при 1-м подходе при передаче значений узлов на границе с р-й зоны на (д + 1)-ю (15) используется дополнительная линейная интерполяция. Кроме того, при 2-м подходе расчет полей в приземной области около здания будет точнее, чем при 1-м (рис. 4). Однако отметим, что поскольку при 1-м подходе сетки перекрываются по одно-
му узлу меньшей зоны, а при 2-м подходе - по двум узлам, то при 2-м подходе число узлов сеток будет несколько больше, чем при 1-м (в одной из задач - на 8,7 %), и, следовательно, время расчета также будет несколько больше.
Таким образом, данные методы оптимизации -метод сгущающихся сеток и метод асинхронного интегрирования по времени - позволяют увеличить точность вычислений в тех местах рассчитываемой области, где она актуальна (вокруг здания), и уменьшить время расчета вследствие уменьшения шагов по расстоянию и по времени в остальной части рассматриваемой области.
д=0 д=1 д=2
1 1
0 |о>- —ь
д=0
2
д=1 д=2
-+ж
ььш
Рис. 4. Наложение сеток в околоземной области: 1-й и 2-й подходы
Примеры расчетов. Была проведена серия расчетов ветровой нагрузки, действующей на различные здания и сооружения, как одиночные, так и в застройках, изложенными выше методами.
Кроме того, в одной из задач было проведено сравнение полученных результатов с результатами натурного эксперимента и с расчетом по формулам СНиП. В данной задаче рассматривалось здание в форме параллелепипеда, размеры которого равны 20*30*25 м, скорость ветра на высоте 100 м (в 4 раза больше высоты здания) равна 6 м/с. Отличие полученных результатов от результатов натурного эксперимента составило 3,4 %, от результатов расчета по формулам СНиП - 4,3 %. Столь малая разница показывает, что изложенные методы достаточно хорошо описывают физические процессы и дают достаточно точные результаты.
Была решена также следующая задача. Рассматривается бесконечная застройка зданий в форме параллелепипеда, в один ряд перпендикулярно направлению ветра. Размеры зданий -16*16*40 м, интервалы между зданиями - 64 м, скорость ветра на высоте 120 м равна 20 м/с. Требуется найти ветровую нагрузку, действующую на здание в данной застройке. Для решения задачи рассматривается область, содержащая одно здание застройки и ограниченная линиями симметрии, проходящими через середины интервалов между зданиями. На линиях симметрии ставятся симметричные или периодические граничные условия. В результате ветровая нагрузка получилась на 12,4 % больше, чем в случае одиночного здания.
Также была решена задача, в которой в рассматриваемой области расположено 2 здания. Они так же, как и в предыдущей задаче, имеют форму параллелепипеда и расположены в ряд перпендикулярно направлению ветра. Размеры зданий и скорость ветра те же, интервал между зданиями - 16 м. Ветровая нагрузка в данной задаче получилась на 20,1 % больше, чем в случае одиночного здания.
Эти два примера обнаруживают существенный недостаток формул СНиП. Ведь такие застройки в СНиП никак не учитываются. А ветровая нагрузка в них значительно возрастает по сравне-
нию с одиночным зданием. И следовательно, требуется повышать устойчивость конструкций этих зданий. Конечно, при расчете конструкций зданий по формулам СНиП их надежность завышается зачастую в разы. Но в таком случае происходит значительное удорожание строительства здания. Однако расчет надежности конструкций сооружений выходит за рамки данного исследования. Целью же данного исследования является разработка методов и программного обеспечения, дающих адекватные результаты, для расчета ветровой нагрузки, действующей на здания и сооружения. Таким образом, разработанный программный комплекс позволяет проводить расчет ветровой нагрузки, действующей на здания и сооружения в различных застройках.
Еще в одной задаче рассматриваются два здания в форме параллелепипеда, которые соприкасаются друг с другом по стороне так, что об-
разуют здание сложной формы. На рис. 5 представлены результаты расчетов данной задачи: стрелками изображено поле скорости. Подобным же образом можно решать задачи на расчет ветровой нагрузки, действующей на здания и сооружения различной формы, представляя их как совокупность нескольких объектов. Заметим, что в СНиП приводится очень ограниченный набор форм зданий и сооружений.
Список литературы
1. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980.
2. Балаев Э.Ф. и др. Численные методы и параллельные вычисления для задач механики жидкости, газа и плазмы: Учеб. пособие. - Иваново, 2003.
3. Строительные нормы и правила 2.01.07-85. Нагрузки и воздействия.
Рис. 5. Результаты расчета ветровой нагрузки, действующей на здание сложной формы: поле скорости Филатов Евгений Юрьевич,
ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», аспирант кафедры программного обеспечения компьютерных систем, телефон (4932) 26-98-29.
