CC BY
145
Ii университета 'ЖУРНАЛ водных / I коммуникации
К. А. Клеванный,
д-р физ.-мат наук,, СПб Центр по гидрометеорологии и мониторингу окружающей среды;
Е. В. Смирнова,
ЗАО «Ленгипроречтранс»
ИСПОЛЬЗОВАНИЕ ПРОГРАММНОГО КОМПЛЕКСА CARDINAL
USING OF MODELING SYSTEM CARDINAL FOR SOLVING HYDRAULIC
PROBLEMS
В статье приводится краткое описание программного комплекса CARDINAL, предназначенного для решения задач гидродинамики поверхностных вод, пример применения комплекса для выбора проектного решения строительства гидротехнических сооружений на р. Вычегда в районе г. Коряжма и Котлаского ЦБК.
A short description of modeling system CARDINAL is given. The system allows to solve different problems offluid dynamics of surface waters. An example of application of the system for design of hydraulic structures on Vychegda river near Koriazhma and Kotlas pulp-and paper mill is given.
Ключевые слова: гидродинамика поверхностных вод, гидродинамическое моделирование, проектирование гидротехнических сооружений
Key words: fluid dynamics of surface waters, hydrodynamic modeling, design of hydraulic structures
m
152]
ИДРОДИНАМИЧЕСКИИ моделирующий комплекс CARDINAL (http://cardinal.x-users.ru) имеет следующие основные функциональные возможности:
• позволяет выполнять расчеты течений, уровня воды, распространения растворенных и взвешенных веществ, температуры и солености, транспорт донных наносов в произвольном водном объекте в двух- и в трехмерном приближениях;
• используются криволинейные гранично-зависимые или прямоугольные координаты и о-координата в вертикальном направлении, обеспечивающая сгущение сетки на мелководье;
• уравнение переноса примеси решается с помощью неявного консервативного конечно-разностного метода высокой точности;
• предусмотрена осушка-затопление ячеек сетки;
• обладает удобным графическим пользовательским интерфейсом;
• пре-процессор, пост-процессор и блок решения объединены в одной программе;
• для совместимости с другими программами используется импорт-экспорт через текстовые файлы;
• возможность импортировать атмосферные данные из метеорологических моделей в формате ОШВ, принятом для передачи метеоинформации;
• возможность задания произвольных аналитических выражений для эмпирических коэффициентов и начальных полей;
• анимация результатов.
Система уравнений, которая решается в этом программном комплексе методом конечных разностей, имеет вид:
ды ды ды ды дс е }др , — + ы — + V — + м— = - е—---— -
■ U--h V--h W-
ôt ôx ôy ôz
ôx p0 Z ôx
1 ÔPa Ô-U Ô-U Ô . ôu (1) ---+ fc V + K (—T + —2-) + — (v ),
ô 2 U ô 2 U ô ôU
Po Ôx
ÔV
Ôx Ôy ÔZ ÔZ '
ÔV
ÔV ÔV ÔV ÔV ÔC g r Ôp ,
— + u — + v — + w — = - g—---— I —— dz -
Ôt Ôx Ôy ÔZ Ôy po Z Ôy
as ^/Ô 2V Ô 2V. Ô , ÔV. (2) a -f„ + rr— +— ) + — (vT—), V-)
Ôx Ôy Ôz ôz
-— P - fc u + K (
po Ôy
3 u д v 3 w
н--+
д х д у
д z
ш
д с dU dV — + — + —= ю ,
dt дх ду s
дТ дТ дТ дТ
--н u--н v--н w — = arsTs
дг дх ду дz
+ K (д2Т , д2Т^ + дТ)
с дх2 ду2 дz с дz
3s 3s 3s 3s
--н u--н v--н w— = msSs +
St дх ду Sz
,32 S д2 S. д , SSn + Kc (—г н—г) н— (v c—)
с \ ^ 2 ^ 2 У -"y \ c S-. s
дх ду 3z 3z
дс дс дс , . дС
--н u--н v--н (w - w0) — =
дг дх ду 3z
. ^ .д2с д2с. д . дс.
= mscs - Яс н кс(—J- + — (v —),
дх ду 3z 3z
(3)
(4)
(5)
(6)
(7)
где u, v, w — декартовы компоненты вектора скорости течения v;
p(T, S) — плотность воды, определяемая по ее температуре T и солености S по зависимости UNESCO;
ро — среднее значение плотности воды (задано 999.7 кг/м3);
g — ускорение силы тяжести; Z — уровень свободной поверхности; PA — атмосферное давление; fc = 2ю sin ф — параметр Кориолиса; ф — географическая широта; ю — угловая скорость вращения Земли; К и vT — коэффициенты горизонтального и вертикального турбулентного обмена соответственно;
U и V — полные потоки; ГО s — объем воды, поступающий в единичный объем за единицу времени из внутренних источников (или изымаемый из него, например, при работе фрезерного землесоса);
Ws — объем воды, поступающей из внутренних источников на единицу поверхности в единицу времени;
с — концентрация растворенных или взвешенных примесей;
w0 — скорость осаждения взвешенных примесей (гидравлическая крупность);
, Т и « — концентрация примеси, температура и соленость воды, поступающей из источников;
X — коэффициент неконсервативно -
сти;
Кс у с — коэффициенты горизонтальной и вертикальной диффузии.
Координата г направлена вертикально
вверх.
На дне при г= -И(х,у) задается касательное напряжение трения Т ь = (т Ь( Х),Т Ь( у)) с помощью квадратичного закона:
д и , _
Ь (X) = Ро^Т д^= "к = ~РоУьи I V |
д V | _
(у) = Ро^Т д^= " * = ~Ро1ъ V 1 V
(8)
где /ь — коэффициент придонного трения.
Для ветровых условий касательное напряжение трения на поверхности Т ^ = (т ж (Х),Т Ж( у)) также задается с помощью квадратичного закона:
3u i _
V„ ^Т ^ = = P„CDW, I W I
3v i —
у) =PoV ^z-=,=PaCDW 1 W |,
(9)
где Са — коэффициент ветрового трения;
^ = (^х ,Жу ) — скорость ветра на высоте 10 м;
ра — плотность воздуха (1.225 кг/м3).
При наличии ледяного покрова касательное напряжение трения на поверхности определяется из соотношений
(10)
X) =~РоАеи I г=с у) = ~Ро1гсе V| V|| г=? > где е — коэффициент трения о лед.
На твердых боковых границах модели нормальная компонента скорости ип равна нулю. Для тангенциальной компоненты ит вводится квадратичный закон трения, аналогичный придонному. Таким образом, граничное условие на твердых боковых границах имеет вид:
00 о-
о
J53|
Un = 0,
K ^T = -fbuAV\ an
(11)
На открытых боковых границах модели граничное условие может быть разных типов. Для первого типа задается нормальная компонента скорости как функция времени: ип = /'ф. Для второго типа задается зависимость уровня воды от времени: С = СО).
На открытых границах, через которые вода поступает в расчетную область, задается временной ход температуры, солености воды и концентрации примеси. На открытых границах, через которые вода вытекает из расчетной области, при условии, что адвекция доминирует над диффузией, граничных условий для с, Т и не требуется.
Для определения коэффициента вертикального турбулентного обмена V возможно использование к-е модели турбулентности, в которой он определяется из соотношения к2
к (12)
v г = с»
£
где к — кинетическая энергия турбулентных пульсаций;
s — скорость диссипации этой энергии за счет внутреннего трения;
с = 0.09.
р
В программном комплексе CARDINAL уравнения для к и s решаются без учета адвективных членов, которые имеют второй порядок малости:
дк д . дк, ЛГЧ ^ = ^ (ут т) + yT (р + N) -е,
dt dz dz
де д ,vT де. е . лгч е
= (— -Г-) + TVT № + Сзе N) - С2— ,
дt дz — dz к к
(13)
где
Р =
f ди \ 2 f а Л ду 2
+
i д z i дz j
N = ■
g Ф,
Prp0 дz
ш
150J
с£ =1.44, =1.92, а =1.3, Pr=0.96.
Для коэффициента с3г при численных экспериментах были подобраны значения, обеспечивающие всплывание вод пониженной плотности и погружение вод с повышенной плотностью: с3^ = 0.001 при устойчивой
2 при неустойчивой.
Для энергии турбулентности и ее диссипации задаются следующие граничные условия на дне (ь) и на поверхности ():
N -о и*2 2 fbub bs = CD PLW 2 Р0
и*3 3/4/ 3/2 Сц К е и* 3/4t3/2 = СМ b3 (14)
еь = KZb Kzb Kzs Kzs
где и* — скорость трения;
к = 0.4 — постоянная Кармана; 2ь и г — параметры шероховатости дна и поверхности соответственно.
Коэффициенты горизонтального турбулентного обмена и диффузии можно определять по формуле Смагоринского:
,, = 0^J if J +(f J - 0,(1 + 0 J.
(15)
Для расчетов эрозии дна и транспорта влекомых наносов к уравнениям системы добавляется уравнение
ht =-
Рр Pg
i
диs + дг1
дх ду
J,
(16)
где pp — плотность частиц в донных отложениях;
pg — плотность донных отложений;
U и V — компоненты вектора расхода донных отложений, которые могут быть определены в CARDINAL из различных эмпирических соотношений.
Сюда входят формулы Мирцхулава-Леви, Гончарова, Эйнштейна, Ван Райна, три вида зависимостей Снищенко-Копалиани, учитывающие высоту донных град, Кнороза, Талмаза и Копалиани для частиц диаметром более 10 мм. Для примера приведем зависимость, полученную Ван Райном:
^ л" f _ V1
U = 0.053,
РР 1 I /3/2 Р - 1 \g d
-1\ D-
стратификации и с3
\\i.P J x V crit , (17)
где d — диаметр частиц грунта;
r _ л
x = PfbU Xcrit =рв P-1 gd, Ds = dt1/3,
,P J
f _ л
ts = :1.44-1012 Рр--1 ?
P ,
0.3
в =
0.24D~\ Ds < 4
0.14 D 0.04D 0.013D 0.055.
-0.64
0.1
0.29
Ds < 10
Ds < 20 D < 150
Эрозия начинается при Т > Т сги . Для расчета скорости осаждения w0 использовалась следующая зависимость, принятая к использованию в ЗАО «Ленгипроречтранс»:
ж = 10.528Ар/2 , (18)
0 18ц + О.б^ТйАр/
Ар = (рр -1000)/1000,ц-
0.01775
-10"
1 + 0.0337í + 0.0002212
где t — температура воды в градусах по Цельсию.
Все значения даны в системе единиц СИ. Сцепление частиц грунта, влияющее на скорость осаждения, не учитывается.
При решении задач в двумерной постановке уровень свободной поверхности воды находится из уравнения (4), а скорости — из решения осредненных по глубине уравнений (1)-(2) (уравнения мелкой воды или уравнения Сен-Венана).
Главное меню программного комплекса состоит из следующих разделов: Project (работа с файлами проектов), Area (создание расчетной области), Depth (ввод поля глубин), Vertical (переход к трехмерной задаче, выбор модели турбулентности, задание числа слоев), Open boundaries (задание открытых границ и условий на них), Atmosphere (задание полей ветра и атмосферного давления, коэффициента ветрового трения), Dynamics (задание начальных полей уровня и скорости, задание коэффициента горизонтального турбулентного обмена), Sources & Pollution (включение расчета распространения примесей, задание источников и начального поля загрязнения), Temperature (аналогично для температуры воды), Salinity (аналогично для солености), Bottom & Ice (описание свойств дна и ледовой поверхности), Devices (установка в расчетной области различного типа «самописцев»), Time (задание шага по времени, расчетного времени и т. п.), Tools (инструменты), Results (вы-
вод на экран результатов расчетов — более 35 видов, возможен просмотр, не останавливая расчет), Window (выбор расположения окон с результатами расчета), RUN (запуск на счет), Help (интерактивная система помощи на английском языке), Print (печать результатов расчета в графический файл или на принтер).
Необходимым этапом разработки гидродинамической модели была ее верификация на аналитических тестах. Сюда входила задача о свободных колебаниях в круговом бассейне, расчет вертикальной структуры ветрового течения в бесконечном канале, расчет течения в канале, вызванном горизонтальным градиентом плотности, расчет переноса примеси в круговом бассейне и др. Все они при дроблении шага сетки дали решения, стремящиеся к точным аналитическим решениям. На рис. 1 показан результат расчета адвекции пятна примеси при нулевой диффузии в круговом бассейне, в котором задано твердотельное вращение жидкости: u = y/R, v= -x/R, где R = 1000 м — радиус бассейна. Начальное пятно задано соотношением c = 100 • exp(-1.5 -10 5 [х - 300)2 + y2}. Как
видно, пятно после двух оборотов возвращается в прежнее положение, практически не изменив своей формы, как это и следует из точного решения. Расчет выполнен в двумерной постановке на сетке 200 х 200 узлов с шагом по времени 60 с.
На рис. 2 показано сравнение результатов расчета течения в канале с резким увеличением ширины с результатами гидравлического моделирования, выполненного во ВНИИГ им. Б. Е. Веденеева. При использовании формулы Смагоринского для коэффициента горизонтального турбулентного обмена (15) без дополнительных настроек модели сразу же была получена длина вихревой зоны, близкая к полученной в лотке, — 5.35 м.
Программный комплекс CARDINAL п
используется для решения разнообразных за- К
дач прикладной гидродинамики и гидравлики. Одна из них связана с выбором оптимального варианта компоновки гидротехнических сооружений. На рис. 3 показан пример расчета поля течений в одном из вариантов проекта строительства нового морского пассажирского терминала в С.-Петербурге.
Рис. 1. Слева — поле скоростей и изолинии концентраций в начальном пятне, справа — изолинии концентраций после двух оборотов пятна в бассейне
Рис. 4. Изолинии скоростей при паводке при наличии мостового перехода, струенаправляющих дамб на обоих берегах и трех пойменных отверстиях в условиях затора в районе проектируемого моста через р. Надым
На рис. 4 показан другой пример — расчетные изолинии скоростей в районе проектируемого моста через р. Надым у г. Надым и г. Старый Надым в русле и на пойме при паводке с расходом воды 14 900 м3/с обеспеченностью 0.33% (1 раз в 300 лет).
Следующим направлением использования программного комплекса CARDINAL
является расчет течений для тренажеров судоводителей фирмой TRANSAS Marine. В частности, тренажер необходим для обучения судоводителей прохождению через су-допропускное отверстие С-1 шириной 200 м комплекса сооружений защиты Санкт-Петер -бурга от наводнений (КЗС). На рис. 5 показан пример поля течений в С-1 при нагонной си-
Рис. 5. Расчет течений в судопропускном сооружении Комплекса сооружений защиты Санкт-Петербурга от наводнений
Рис. 6. Вертикальный разрез концентраций взвешенных веществ при добыче железно-марганцевых конкреций в карьере «Кургальский» через 12 часов работы землесоса. Западный ветер 0.5 м/с
туации. Расчет сделан на прямоугольной сетке с шагом 6 м.
Большое количество различных расчетов выполнялось для оценки экологического воздействия при различного типа гидротехнических работах и при прогнозе последствий от аварийных разливов и выбросов. Результаты этих расчетов (объемы воды, в которых концентрации взвешенных веществ превышают заданные значения и площади дна, на ко-
торых образуется слой осадка толще заданной величины) используются для расчета ущерба ихтиофауне. На рис. 6 — вертикальный разрез концентраций взвешенных частиц на карьере «Кургальский» в Финском заливе, к северу от о. Мощный, сбрасываемых с судна при добыче железно-марганцевых конкреций. Судно ходит по участку галсами вперед и назад, поднимая грунт с конкрециями, а пустую породу сбрасывает в залив. На рис. 6 видны
ш
158]
Рис. 7. Концентрации Сб-137 на дне Копорской губы через 10 суток после выброса в атмосферу при северо-восточном ветре
Рис. 8. Сравнение прогностического и фактического хода уровня на в/п Горный институт
за период 26.10-7.12.2006 г.
поля мутности, образованные при настоящем и предыдущем галсах. Расчеты выполнялись раздельно для каждой из фракций, затем результаты складывались. Видно более быстрое выпадение тяжелой фракции.
Пример расчета оседания на дно Копор-ской губы радиоактивного изотопа Cs-137 после аварийного выброса в атмосферу с ЛАЭС-2 приведен на рис. 7. В данном случае в CARDINAL импортировались данные о концентрациях Cs-137 на поверхности из атмосферной программы ZONA. Учитывалась адсорбция изотопа на взвешенных частицах. На рисунке результаты расчета при северо-восточном ветре 1.5 м/с.
Пожалуй, наиболее важным является использование программного комплекса CARDINAL для оперативного прогноза
уровня воды в Балтийском море и, в частности, в Санкт-Петербурге. Система прогноза полностью автоматизирована и работает в СЗУГМС с декабря 1999 г. Прогноз метеопараметров над Балтийским морем поступает ежедневно (а с апреля 2008 г. — четыре раза в сутки) из модели ШЯЬАМ Шведского гидрометеорологического института. Для прогноза уровня используется детальная модель Балтийского моря. Заблаговременность прогноза — 48 часов. Пример сопоставления расчетного и фактического хода уровня воды в С.-Петербурге представлен на рис. 8. Соответствие достаточно хорошее. Постоянно ведется работа по повышению точности прогноза. После завершения строительства КЗС прогноз будет необходим для управления затворами КЗС.
Список литературы
1. Аверкиев А. С., Клеванный К. А. Определение траекторий и скоростей циклонов, приводя- ы щих к максимальным подъемам воды в Финском заливе // Метеорология и гидрология. — № 8. К
— 2007. С. 55-63. 4ГГ1
2. Вольцингер Н. E., Клеванный К. A., Пелиновский E. Н. Длинноволновая динамика при-^^59 брежной зоны. — Л.: Гидрометеоиздат, 1989.
3. Клеванный К. А., Матвеев Г. В. Моделирование воздействия гидротехнических сооружений на гидрологический режим восточной части Финского залива с помощью программной системы CARDINAL // Финский залив-96: материалы Международного научно-практического симпозиума. — СПб.: Изд-во «Балтвод», 1996. — С. 51-58.
ш
160J
4. Клеванный К. А. Наводнения в С.-Петербурге при незавершенном комплексе защитных сооружений // Известия АН. Сер. географическая. — № 2. — 2002. С. 80-88.
5. Клеванный К. А. Прогноз изменений в режиме течений и концентрациях загрязняющих веществ в Невской губе, связанных с окончанием строительства сооружений защиты С.-Петербурга от наводнений // Акватерра: тр. V Международной конференции, 12-15 ноября 2002 г. — СПб.: Водоканал, 2002. — С. 81-83.
6. Прогностическая оценка дисперсии радионуклидов в крупных водоемах на примере озера Эворон // ЭРАС-3. Экология регионов атомных станций: сб. ст. / A. A. Апухтин [и др.]; под ред. Ю. A. Егорова. — М.: Изд-во Атомэнергопроект, 1995. — Вып. 3. С. 134-145.
7. Расчет гидрологического режима Невской губы / Н. Е. Вольцингер [и др.] // Метеорология и гидрология. — 1990. — № 1. С. 70-77.
8. Alkhimenko A., Bolshev A., Jalovlev A., Klevannyy K., Liukkonen S. Modeling oil pollution under ice cover // Proc. Vll-th International Offshore and Polar Engineering Conference. — Honolulu, USA: International Society of Offshore and Polar Engineers, 1997. — Р. 594-601.
9. Klevannyy K. A., Preobrazhensky L. Yu., Voltzinger N. E., Zolnikov A. V. Body-fitted co-ordinates for calculating long-wave disturbances and impacts of coastal engineering structures on the hydrodynamic and pollution regimes // Proc. of the Ill-th International Conference on Coastal & Port Engineering in Developing Countries. — Mombasa, Kenya, 1991. — Vol. 2. — Р. 1311-1320.
10. Klevannyy K. A., Matveyev G. V., Voltzinger N. E. An integrated modelling system for coastal area dynamics // Int. J. for Numerical Methods in Fluids. — 1994. — Vol. 19. — Р. 181-206.
11. Klevannyy K. A. Simulation of storm surges in the Baltic Sea using an integrated modelling system CARDINAL // Proc. of the XIX-th Conference of the Baltic oceanographers. — Sopot, Poland, 1994. — Р. 328-336.
12. Klevannyy K. A., Matveyev G. V. Integrated modeling system CARDINAL and its applications to estimation of environmental impacts of coastal engineering in the Gulf of Finland // Proc. of the IV-th International Conference on Coastal & Port Engineering in Developing Countries. — Rio-De-Janeiro, Brazil, 1995. — Vol. 3. — Р. 2120-2134.
13. Klevannyy K. A. Modelling floods for the present state of St. Petersburg Flood Protection Barrier // Proc. of Second Study Conference on BALTEX. — Juliusruh, Island of Rogen, Germany, 1998. — Р. 106-108.
14. Klevannyy K. A. Protection of St. Petersburg against floods with uncompleted barrier: mathematical model study // Proc. of the V-th International Conference on Coastal & Port Engineering in Developing Countries. — Cape Town, South Africa, 1999. — Р. 626-637.
15. Klevannyy K. A., Gubareva V. P., Mostamandi M. S. W., Ozerova L. B. Water level forecasts for the Eastern Gulf of Finland // Bulletin of the Maritime Institute in Gdansk. — 2001. — Vol. 28. — № 2.
— Р. 71-87.
16. Klevannyy K. A., Smirnova E. V. Simulation of current and water pollution changes in the Neva Bay after completion of St. Petersburg flood Protection Barrier // Environmental and chemical physics.
— 2002. — Vol. 24. — № 3. — Р. 144-150.
17. Klevannyy K. A. Simulation of impact of flood protection barrier on Neva Bay water regime // Sixth International Conference on Coastal & Port Engineering in Developing Countries. Colombo, Sri Lanka, Book of Abstracts, 2003. pp.137 - 138.
18. Klevannyy K. A., Gliantseva O. V., Mostamandi M. S. W. Evaluation of new formulation of wind drag in the operational model of the Baltic Sea // Proceedings of the 6th HIROMB Scientific Workshop, St. Petersbvurg, 8-10 Sept. 2003. — St. Petersbvurg, 2003. — Р. 21-32.
