автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Моделирование гидродинамических процессов в мелководных водоемах на оптимальных криволинейных сетках

На правах рукописи

сГ

Гч.

ВАСИЛЬЕВ ВЛАДИСЛАВ СЕРГЕЕВИЧ

МОДЕЛИРОВАНИЕ ГИДРОДИНАМИЧЕСКИХ ПРОЦЕССОВ В МЕЛКОВОДНЫХ ВОДОЕМАХ НА ОПТИМАЛЬНЫХ КРИВОЛИНЕЙНЫХ СЕТКАХ

05.13.16- применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Ростов-на-Дону 1997

Работа выполнена в Таганрогском государственном радиотехническом университете.

Научные руководители: действительный член АЕН РФ,

доктор технических наук, профессор Л.С.Берштейн доктор физико-математических наук,

профессор А.И.Сухинов

Официальные оппоненты:

доктор физико-математических наук,

профессор Л.А.Крукиер доктор физико-математических наук,

профессор А.А.Илюхин

Ведущая организация:

Институт математического моделирования РАН, Москва

Защита состоится "<50 " окюя&з 1997 г. в И часов на заседании диссертационного совета K063.52.I2 по физико-математическим и техническим наукам в Ростовском государственном университете по адресу: 344090, г.Ростов-на-Дону, пр.Стачки, 200/1, корпус 2, вычислительный центр РГУ.

С диссертацией можно ознакомиться в научной библиотеке РГУ по адресу: ул.Пушкинская, 148.

Автореферат разослан "2Э" сентября 1997 г.

Ученый секретарь диссертационного совета,

кандидат физико-математических наук Г.В.Муратова

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность тепы диссертации. Большинство внутренних водоемов России и многие прибрежные районы морей являются мелководными водоемами, для которых отношение характерных размеров в горизонтальных направлениях составляет 104 и более. Но именно эти водные экосистемы испытывают наибольшую антроттогрнт'к> нагрузку, связанную со строительством портов, прокладкой судоходных каналов, а также с различного рода загрязнениями. Многие из этих экосистем являются уникальными и весьма "хрупкими" (такие как Таганрогский залив и Азовское море в целом, северный Каспий и др.).

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

Среди задач водной экологии первичными являются расчеты гидрофизических параметров. Гидродинамические характеристики водоемов (поля скоростей водного потока, функции возвышения уровня) являются входными данными для расчета распространения загрязняющих веществ и развития всего живого в водной среде, начиная от фитопланктона и заканчивая млекопитающими.

Реальные водоемы представляют собой области сложной формы с объективно предопределенной границей. Граничные узлы прямоугольной сетки декартовых координат, покрывающей водоем, могут точно не попадать на контур береговой линии акватории. Это приводит к тому, что краевые условия ставятся в точках, смещенных от реальной границы, либо, при усечении граничных ячеек сетки, вблизи границы концентрируются неоднородности, связанные с неравномерностью шагов по пространственным направлениям. Кроме того, дно мелководных водоемов может быть изрезано углублениями кусочно-линейной или криволинейной формы, не совпадающей с прямоугольными декартовыми координатными линиями, естественного и искусственного происхождения (бывшие русла рек, прорытые судоходные каналы и т.д.), наличие которых может во многом определять картину течений: отдельные линии тока могут повторять форму этих углублений. Сравнение аппроксимаций различных кривых показывает, что для адекватного приближения их формы ломаными, вершины которых лежат на кривых, требуется в 5...10 раз меньше отрезков по сравнению со ступенчатыми ломанными, смежные стороны которых взаимно перпендикулярны или сонаправлены. Поэтому

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

Цели и задачи исследования. I. Построение гидродинамической модели, удовлетворящей основным законам сохранения (массы, импульса, энергии) и учитывающей испарение воды и выпадение осадков, поскольку для объектов математического моделирования -Азовского моря и Таганрогского залива объем испаряющейся жидкости (преимущественно в летние месяцы) сравним с годовым речным стоком, а годовой объем выпадающих осадков уравновешивает водообмен между Азовским и Черным морями.

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

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

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

Интегрированием по вертикальной координате трехмерных уравнений Навье-Стокса и неразрывности для вязкой несжимаемой жидкости в однородном поле тяжести во вращающейся системе отсчета получено семейство двумерных гидродинамических моделей мелководного водоема, учитывающих испарение еоды и выпадение осадков. Следствием уравнений баланса импульса моделей являются уравнения баланса полной механической энергии, принимающие дивергентную форму в силу уравнений баланса массы (энергетическая нейтральность) и обладающие строгой диссипативностью за счет действия сил внутреннего вязкого трения, в частности, над поверхностью (дна), описываемой гармонической функцией [1-3].

Многие уравнения математической физики содержат оператор Лапласа, форма которого с точностью до множителя инвариантна относительно конформных преобразований координат, минимизирующих (в случае их существования) значение функционала Дирихле в

рассматриваемой области. Это широко используется в методах построения как регулярных, так и нерегулярных сеток. В предшествующих данному исследованию работах была продемонстрирована высокая эффективность методов построения сеток для областей сложной формы, основанных на минимизации в расчетной области непосредственно значения функционала качества требуемой сетки, обеспечивающей независимость от выбора начального приближения. Однако, из-за чрезвычайной сложности точного решения задачи линейного поиска для функционала Дирихле в расчетной области в методах 1-го порядка (К.Сагса111е1, Б.Е.Кеппоп, С.З.Ш11кга71с1г) не использовался функционал Дирихле. А в имеющихся методах 2-го порядка (С.А.Иваненко) элементы матрицы Гессе задачи оптимизации менялись скачкообразно, что сказывалось на работе алгоритма.

В настоящей работе за счет усовершенствований, устраняющих скачкообразный характер изменения элементов матрицы Гессе, удалось повысить эффективность имеющихся методов 2-го порядка [4]. Благодаря замене при точном решении задачи линейного поиска исходного функционала имеющим решение приближением, гладко передающим (до производных 2-го порядка включительно) поведение исходного функционала, удалось разработать эффективные методы 1-го порядка, направленные на получение невырожденного преобразования координат 15). Используя неравенства Адамара и Коши, удалось расширить и обобщить классы используемых функционалов [4,53.

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

Для оценки решений параболических и эллиптических уравнений, главным образом, в системах координат, имеющих особенности, доказана теорема сравнения, оценивающая норму решения в равномерной метрике С выражением, содержащим норму С правой части и дискретный аналог интеграла от метрических коэффициентов монотонного сеточного оператора.

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

скорости сходимости при определенном расщеплении оператора (определенной локально-двумерной схеме (ЛДС)) являются неу-лучшаемыми. При доказательстве не использовалось разделение переменных и построение мажорирующих решений в явном виде (что может быть затруднительно или невозмохно) [6,7].

Практическая значимость. Разработанные метода построения сеток могут быть использованы для любых приложений, требующих отображения заданной области в криволинейную область более удобной формы. Построенные в результате оптимизации функционала Дирихле отображения на акваторию Азовского моря и Таганрогского залива Г-образной области (море - 50*20, залив - 27'10 узлов) и прямоугольников 83*20 и 80*20 узлов максимально точно передают очертания береговой линии. По полученным на основе многолетних наблюдений 812 отметкам глубин форма дна оптимальным образом приближена интерполяционно-сглаживающей поверхностью, близкой к гармонической сеточной функции, что для песчано-илистого дна Азовского моря и Таганрогского залива является приемлемым. И сетки, и поверхность дна могут служить составной частью дискретных моделей других процессов в Азовском море и Таганрогском заливе. Реализованный алгоритм численного расчета ветровых течений с высокой точностью удовлетворяет закону сохранения массы [1-3,81. Все программы доведены до работоспособного состояния.

Достоверность полученных результатов основана на строгих математических выводах (при отмечаемых ограничениях) всех положений и максимально точном учете имеющихся входных данных, а также па качоитишшш, а в ряде случаои и количественном соответствии результатов моделирования реально наблюдаемым течениям в Азовском море и Таганрогском заливе.

Математический аппарат. Полученные результаты основаны на использовании законов классической механики, концепции сплошной среды, уравнений математической физики, в том числе уравнений Навье-Стокса и неразрывности, дифференциального и интегрального исчисления функций многих переменных, матричной алгебры, тензорного исчисления, вариационного исчисления, методов оптимизации (нелинейного программирования), методов расщепления, метода конечных элементов, численных методов (решения систем линейных алгебраических уравнений, определения действительных корней полиномов), теории устойчивости разностных уравнений А.А.Самарского.

Аппробация. Результаты работы докладывались на - XLIII Всесоюзной научной сессии, посвященной Дню радио

(Москва, 1988),

- Шкоде-семинаре по комплексам программ математической физуши (Ростов-на-Дону, 1990),

- IV Всероссийской школе молодых ученых "Численные метода механики сплошной среды" (Абрау-Дюрсо, 1992),

- IV Всесоюзном совещании "Проблемы построения сеток для решения задач математической физики" (Сунгуль (Челябинск-70), 1992),

- Областной научно-технической конференции, посвященной Дню радио (Ростов-на-Дону, 1993),

- Международной научной конференции, посвященной 100-летию со дня рождения Н.Г.Чеботарева (Казань, 1994),

- Всероссийском симпозиуме "Математическое моделирование и компьютерные технологии" (Кисловодск, 1997 ),

- Международной конференции "Математические модели физических процессов и их свойства" (Таганрог, 1997),

а также на научных семинарах кафедр вычислительной математики и вычислительного эксперимента, информатики, на конференциях профессорско-преподавательского состава Таганрогского государственного радиотехнического университета и на семинарах лаборатории вычислительного эксперимента вычислительного центра Ростовского государственного университета.

Публикации. По материалам диссертации опубликовано девять печатных работ.

Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения и списка использованных источников, изложена на 147 страницах машинописного текста, включает 40 рисунков и 222 библиографических названия.

С0ДЕР1АНИЕ РАБОТЫ

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

Первая глава посвящена гидродинамическим моделям мелководного водоема, состоит из 6-и параграфов и предваряется обзором имеющихся моделей и методов решения уравнений Навье-Стокса.

В §1 интегрированием уравнений Навье-Стокса и неразрывности для вязкой (в линейном приближении) несжимаемой (плотность Р^сопэг) жидкости во вращающейся с постоянной угловой скоростью

Q, составляющей угол -ö с вертикалью, системе отсчета в однородном поле ускорения свободного падения g по вертикальной координате 2 от -h до С получены уравнения гидродинамической модели, учитывающей испарение жидкости и выпадение осадков ((и/р) - испаряющийся - выпадающий в виде осадков в единицу времени слой)

VyE*Vfta/p;=0' (1)

и'Асиии2/н) > К л/«) /

+ (тYp) [ш-CJU/H] AC-ubA/ij +F*gx+F*bx+2ífísln-ü, (2)

VW > [cuJ2/H] 'y+^^v РУЯ)

Hr/p) [aV-Gv (V/H) ЛС-vbLh] +ry+F*bv-2Wsiirt, (3)

где слагаемые, имеющие вид и размерность граничных вязких напряжений на поверхности жидкости отнесены на счет обобщенной силы трения ветра о поверхность F*, а на дне - на счет обобщенной силы трения о дно Ff,

£ I

U= ¡udz, V= |vdz, иии- горизонтальные компоненты вектора -к -н

скорости жидкости в точке (x,y,z) в момент времени f, h - высота столба жидкости под невозмущенной поверхностью, С - возмущение свободной поверхности, H=h+С - полная глубина, Cuu' ^uv' ^vv' Cu* ^v ~ упрзвля10^6 коэффициенты:

С С С

fu2dz=C H'1U2, Гuvdz=C H'1UV, Гv2üz=C Н~1Ч2,

¿ uu i uv ¿ vv

-h - h -h

U -С И'111, v -П И' V, П MC >1, С2 С ,

eu а V uu' uu uu uu w

t¡ - первый коэффициент вязкости, Л=дг/дх2+д2/ду2 - двумерный оператор Лапласа.

В §2 показано, что для упрощенной (CuusCuv)=Cvij=7) модели (I)-(3) над неподвижным h't=0 дном выполняется закон сохранения полной механической энергии - суммы потенциальной энергии ü=gH(4~h)/2 в результирующем поле тяжести и положительно определенной квадратичной формы К=Е- [иг^г)/(2Н) интегралов U и У (аналога кинетической энергии) в дифференциальной форме, где E=E(x,y,t)>0 удовлетворяет уравнению переноса (в частности Est)

E't + (U/H )E'x+(V/U }Ey=0.

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

Модель оказывается строго диссипативной за счет действия сил внутреннего вязкого трения над поверхностью (дна), являющей-

ся гармонической функцией

АЛ=0. (4)

В §3 выводятся уравнения других гидродинамических моделей, допускающих существование уравнений баланса полной механической энергии. В частности показано, ценой каких упрощающих предположений, внедряемых на более ранней стадии вывода модели (I)-(3), получается модель с диссипативным фильтром, встречающимся у других авторов (Н.Е. Вольцингер и др.).

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

В §4 проводится дискретизация модели (1)-(3). В частности, обосновывается выбор естественных (скорость - возвышение уровня) переменных, криволинейных сеток, конечно-элементной пространственной' яттроксгмяцт* И М°ТПДЯ ря^отт.лот'а ттп фи'эинагчп«" ттпп-

цессам - аналога метода поправки к давлению при дискретизации по времени, а также отнесение на новый или на предыдущий временной слои различных слагаемых уравнений (1)-(3). После отбрасывания (ввиду малости вклада) слагаемых, содержащих вторые производные от Ь, Я и С. базисные функции выбираются из условий минимально необходимой гладкости. При этом различные сеточные функции оказываются заданными на различных компонентах (узлах, ребрах, ячейках) сетки (разнесенная сетка). В результате, на т - шаг по времени или итерационный параметр накладывается ограничение типа критерия Куранта.

Задача решается в интегралах У и У (а не в усредненных по полной глубине компонентах) от декартовых (а не от контра-вариантных) компонент вектора скорости.

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

деленный оператор при неотрицательной замороженной полной глубине H предыдущей итерации.

В §5 приводится конечно-элементная пространственная аппроксимация для отдельных слагаемых схемы.

В §6 излагается постановка условий на границе и неявный итерационный метод решения задачи для возвышения уровня.

В граничных узлах сетки задаются значения функций U, V и О, V (вязкая жидкость - условия прилипания - задача Дирихле), что эквивалентно заданию на границе значений производных Ç^ и для возвышения уровня Ç (задача Неймана).

Для нахождения решения ç на новом слое используется итерационный процесс [91, а на каждой итерации по нелинейности получающаяся система линейных алгебраических уравнений решается методом Гаусса-Зейделя, поскольку структура блочно-трехдоа-гональной матрицы оператора такова, что в диагональных блоках на оптимальной сетке имеет место значительное диагональное преобладание (на идеальной сетке из квадратных ячеек диагональные блоки будут диагональными матрицами), и, кроме того, на основе метода Гаусса-Зейделя предложен и обоснован (С.С.Маханов, А.Ю.Семенов) численный алгоритм сквозного счета для моделирования двумерных течений жидкости, обеспечивающий строгую неотрицательность получаемых глубин для любых режимов течения.

Решения U, V, а после определения С. и Û, V находятся по явным схемам.

Вторая глава посвящена методам построения невырожденных криволинейных расчетных сеток для областей сложной формы, состоит из 4-х параграфов (§7-§10) и предваряется обзором имеющихся методов построения регулярных и нерегулярных, в т.ч. связанных с решением и подвижных сеток.

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

В §8 на основании неравенств Адамара и Коши получены неравенства (парные греческие индексы - суммирование от 1 до п)

о < 1-тег\к 1-ва1+ KJ-gat| « ^ ¡¡-z ilgai, Ktjm.

" а а " i=1 а 1=1 а

0 < K«.gai| S S nftK'Sar W-J*1'

■ 1 i = I t — I

где К*., Ki.J^n - компоненты смешанного тензора 2-го ранга,

и , , - компоненты метрического и ассоциированного метрического тензоров преобразования координат

(х\х?,... ,хп) н-у (5)

связывающие определители и следы с произведениями диагональных элементов положительно определенных матриц симметризованных тензоров 2-го ранга, равных сумме тензоров, ассоциированных смешанному тензору К, описывающему свойства анизотропной среды, и тензору Кт, полученному транспонированием К, и на основании их предложены семейства функционалов качества преобразований (5)

<3£{ ас£

{= I

Я---Р

= XX-Я-.-ТР = Я---Р

п а^ ас _ _ . .

2 КА—»--Л- п кй--»--к+кг|

Л дха дх? 1=1 Р дха вхР I

сЬг'йх2.. .(±гп =

тг ™, £ .

{=; к

, п к* 1=1 а

к+кт\/8

(6)

п с ае7 ае5 п с, ас7 б

у к0--^1----п ка—

¿1 е'1 дха Г' • р

¡к+кг|/в-

дх$ 1=1 * (Зх'йх2...^.

дха

аг

дхр

(7)

Отмечено, что в случае изотропной среды и п=2;3 семейства содержат многие широко используемые функционалы (Дирихле, А.Ф.Сидорова, Г.П.Прокопова, методов Врэкбилла-Зальцмена и "упругой паутины", В.Д.Лисейкина и др.).

Показано, что в случаях Р(и^,и>)=Р(у/,ю) (а если и

Не Зависят 0,г то и Р(и,У,и>)=Р(И/и>)) минимумы

функционалов (6),(7) (в случае существования) достигаются на преобразованиях координат (5), удовлетворяющих обобщенным условиям ортогональности (при совпадающими с обычными)

.. сЯ{ д^ , .. дха

=[к^з-—■— = о- ни.*.

а г ас

"дх^'дх^ с р а) д^ д^

При этом, в случае функционала (6) эллиптический дифференциальный оператор 2-го порядка

СИУСК^айФ;

д

дх*

еф"

ь-а__

Р а^Р

зф а

д^ дха

к!

дх?

Рдх? дха

д^Ф

не содержит смешанных производных (а в случае однородной

ак^/ах1=о, ак^/асас 1/агт^о,

среды при Р(и,ь,ю)=и и К^ , - и кососимметрической со-

ставляющей, что можно рассматривать как обобщение конформных

Рис.2. Изолинии поверхности дна ( &.П = 1 м )

Рис.3. Поле течений (западный ветер)

Рис.4. Ветровые нагоны (западный ветер, аС = ? см)

преобразований координат).

С другой стороны, при P(u,v,w)=P(w) минимумы функционалов (6) и (7) достигаются на обобщениях

JК+КтJ/g = const (при Р(для функционала (6)),

gJК+К2*| = const (при P(u,v,w)/A/Yw+B для функционала (7)), условий эквиареальности (J=const) преобразования координат (Б),

где g=det|gtJJ, jK^KT|=det|KjfK^|, <4=const, B=const.

Для случаев однородной и изотропной среды приводятся уравнения Эйлера-Лагранжа, сложность которых свидетельствует о том, что методы, основанные на минимизации непосредственно значения функционала, могут не уступать в эффективности методам, основанным на решении систем алгебраических уравнений).

В §9 среди семейств (6) и (7) выделяется подмножество

I =

(8> j

s'y' - х'и',

содержащее как наиболее важных представителей - функционалы

= n[su+i?2}<mi = Jp™^

"s22. . [ÇSug

(10)

1 - ft^*» - iPf^.

I = JJg"'ctKjy = (II)

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

Функционалы (9)-(П) объединяет наличие якобиана J в знаменателях подынтегральных функций, служащее барьером при переходе в область вырождения (J^O), что монет быть использовано в методах построения невырожденных сеток.

В работах С.А.Иваненко предложена замена J на

J*= max(J,e), s>0, (12)

однако трудности, встречавшиеся при построении сеток, заставили искать пути совершенствования этого подхода.

Рассмотрение замены (12) с позиций функций штрафа означает, что величина штрафа не зависит от положения решения в области вырождения (<7<0, точнее J<s, е>0). Рассматривая функцию f(J)~7/е как многочлен Тейлора 0-го порядка для функции f(J)=l/J в точке J=e, предлагается заменить барьерную функцию f(J)=1/J при J^e функцией штрафа, где величина штрафа определяется значением многочлена Тейлора для функции f(J)=l/J в точке J=e более

- 15 -

высокого порядка (1-го, 2-го и т.д.)

к=0

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

В §10 построен метод наискорейшего спуска для подмножества (8). Наибольшие трудности при этом вызывает точное решение задачи линейного поиска на каждом шаге. Предлагается за направление поиска выбирать направление градиента дискретного аналога

(13)

а

функционала из подмножества (8),

где <р[г= зЦд^.Л^.Л^АГ.]./^), Jlr= МДГГ ЛХ^,

'Ы = Н'^г/^Н^хЛ)2)/^Г' ^[./^.е)

(суммирование ведется по всем четырем углам (а) всех ячеек сетки, а индексами I и г отмечены величины, относящиеся к "левому" и "правому" ребрам каждого угла), а за величину шага - положение минимума специальным образом формируемого полинома 4-й степени

где ф1гСО = +

Jгr(t■) = [л^г Д^д] [ЛУ1+ ЛуГ^] - [дХ1+ ЛУГД] ,

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

Алгоритм обобщен на случай методов сопряженных направлений. Демонстрируются достигнутые эффекты.

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

В §11 приводится построенная градиентными методами §10 сетка для Азовского моря и Таганрогского залива (море - 50*20, залив - 27«10 узлов) (рис.1).

X и —

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

1®= юг® + П-чоД® (14)

помимо взятой с весом и (ОФК1) меры интерполяции I® - суммы квадратов разностей между измеренными и интерполяционными значениями глубин присутствует слагаемое I® - мера гладкости -дискретный аналог функционала Дирихле

V = б22И2]ему =

= Я(б22И2- бмК)2]^-

необходимое условие минимума которого представляет собой соотношение (4), обеспечивающее строгую диссипативность моделей §1 и §3 за счет действия сил внутреннего вязкого трения.

Функционал (14) минимизируется методом наискорейшего спуска. На каждом шаге задача линейного поиска допускает точное решение.

При выборе веса интерполяции о>$0.96 на полученных поверхностях дна отсутствуют 13-метровые глубины, а при выборе ш>0.98 получаемая поверхность не адекватна действительности. Оптимальное значение находится в пределах 0.965({^0.9Ч0. При таком выборе веса меры интерполяции время сходимости процесса минимизации функционала (14) не превосходит или сравнимо со временем расчета метрических коэффициентов и локальных координат отметок измерения глубин.

В §12 производится переход к безразмерным переменным и приводятся результаты моделирования ветровых течений при различных скоростях и направлениях ветра.

Четвертая глава посвящена методам получения оценок решений сеточных уравнений и состоит из 5-и параграфов (§13-§17).

В §13 доказывается играющая важную роль в дальнейшем

теорема сравнения. Для задачи

Ах>к.1=^ь,1' ^О'К-1- 1=1 икг=0, 1=1,Ь,

Г х К. Гьъ~1. г К. ». г)Кг+

+ К 1-1/2К г-т]+аь. г* 1/гК Г"*, г

к=1,К-1, 1=1,Ъ, Лио. ГЬ1/2. г К. Г"». 1)ПТ2+[ао. 1-1/2^0, г "ио. 1-т) + ^о.т/гко.^о.т})*?'

К г/г. г/ьъ-иг. г}^{ъъ+1/2. г}/г) •

где А - положительная константа, справедлива оценка |ы| где (1У|с=тах^шах|ш^

В §14 дается постановка начально-краевой задачи для уравнения теплопроводности в тороидальных координатах.

В $15 для поставленной задачи строится локально-двумерная схема.

В §16 на основе доказанной (§13) теоремы устанавливается сходимость решения разностной (§15) задачи к решению начально-

краевой (§14) задачи при т-ьО, со скоростью

при а=0 становящаяся неулучшаемой.

В §17 анализируется скорость сходимости других локально-двумерных схем и проводится сравнение с локально-одномерными.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ, ПРЕДСТАВЛЕННЫЕ К ЗАЩИТЕ

1. Построено семейство двумерных гидродинамических моделей мелководного водоема, учитывающих испарение жидкости и выпадение осадков.

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

2. Предложены семейства функционалов качества преобразований координат. Показано, что семейства содержат функционалы, минимумы которых (в случае существования) достигаются на обобще-

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

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

3. Для акваторий Азовского моря и Таганрогского залива построены оптимальные криволинейные сетки. Форма дна приближена интерполяционно-сглаживающей поверхностью, близкой к гармонической функции (рис.2). Получены поля течений (рис.3) и нагонов (рис.4) при различных направлениях и скоростях ветра.

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

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

Основные результаты диссертации опубликованы в следующих работах

1. Сухинов А.И., Васильев B.C. Гидродинамическая модель мелководных протяженных водоемов и ее численная реализация /Алгебра и анализ. Тезисы докладов Международной научной конференции, посвященной 100-летию со дня рождения Н.Г.Чеботарева.- Казань: Издательство Казанского университета, 1994, Ч.II, С.127-128.

2. Сухинов А.И., Васильев B.C. Расчет ветровых течений в Азовском море на основе новой модели "мелкой воды" методом конечных элементов на оптимальных сетках /"Математическое моделирование эколого-экономических систем". Тезисы докладов на Всероссийском симпозиуме "Математическое моделирование и компьютерные технологии".24-26 агагреля 1997 г.-Кисловодск: Издательство Кисловодского института экономики и права, 1997, С.95-97.

3. Васильев B.C. Двумерные модели гидродинамики мелководного водоема на криволинейных сетках из выпуклых четырехугольников/ "Математические модели физических процессов и их свойства". Тезисы докладов Международной конференции (30 мая - 3 июня 1997 года). - Таганрог: Таганрогский государственный педагогический ин'"1 и-!j•!, 1007. - С.94.

4. Васильев B.C., Фрадшш Ь.Г. Численный алгоритм задачи

расчета сеток для параллельных ВС со структурной реализацией вычислений /Численные методы механики сплошной среды. Тезисы докладов IV Всероссийской школы молодых ученых. - Красноярск, 1992, С.19-21.

5. Сухинов А.И., Васильев B.C. Применение алгоритмов построения сеток для областей сложной формы /Численные методы механики сплошной среды. Тезисы докладов" IV Всероссийской школы молодых ученых. - Красноярск, 1992, С.100-102.

6. Сухинов А.И., Васильев B.C. Локально-двумерные схемы для решения трехмерного уравнения теплопроводности в тороидальных координатах на многопроцессорных системах /XLIII Всесоюзная научная сессия, посвященная Дню радио. Тезисы докладов, - М.: Радио и связь, 1988, С.78.

7. Сухинов А.И., Васильев B.C. Локально-двумерные схемы для аппроксимации трехмерного уравнения теплопроводности в тороидальных координатах //Изв. вузов. Математика, - 1996, J63, С.58-67.

8. Сухинов А.И., Васильев B.C. О моделировании ветровых течений в Таганрогском заливе с использованием квазиоптимальных сеток /Областная научно-техническая конференция, посвященная Дню радио. - Ростов-на-Дону, 1993, С.26.

9. Сухинов А.И., Васильев B.C. Математическое моделирование переноса зарядов в приборах с зарядовой связью //Дифференц. ур-ния, 1992, Т.28, Х2, 0.346-35"

/

¿V' „ 4 / ?

ТАГАНРОГСКИЙ ГОСУДАРСТВЕННЫЙ РАДИОТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

На правах рукописи

Васильев Владислав Сергеевич

У

МОДЕЛИРОВАНИЕ ГИДРОДИНАМИЧЕСКИХ ПРОЦЕССОВ В МЕЛКОВОДНЫХ ВОДОЕМАХ НА ОПТИМАЛЬНЫХ КРИВОЛИНЕЙНЫХ СЕТКАХ

05.13.16 - Применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

Диссертация на соискание ученой степени кандидата физико-математических наук

Научные руководители: академик АЕН РФ, д.т.н., проф. ^/^у^у Л.С.Берштейн,

Таганрог, 1997 г.

- 2 -ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ.......................................................3

ГЛАВА I. ГИДРОДИНАМИЧЕСКИЕ МОДЕЛИ МЕЛКОВОДНОГО ВОДОЕМА........10

§1. Вывод уравнений гидродинамической модели................II

§2. Уравнение баланса аналога полной механической энергии..Л8

§3. Другие гидродинамические модели мелководного водоема----22

§4. Дискретизация гидродинамической модели..................24

§5. Конечно-элементная аппроксимация слагаемых схемы........30

§6. Постановка граничных условий и решение уравнения для

возвышения уровня........................................39

ГЛАВА II. ПОСТРОЕНИЕ КРИВОЛИНЕЙНЫХ СЕТОК ИЗ ВЫПУКЛЫХ

ЧЕТЫРЕХУГОЛЬНИКОВ........................................46

§7. Вспомогательные соотношения.............................47

§8. Семейства функционалов..................................52

§9. Квазиньютоновские методы................................65

§10. Градиентные методы наискорейшего спуска и сопряженных

направлений..............................................74

ГЛАВА III. МОДЕЛИРОВАНИЕ ВЕТРОВЫХ ТЕЧЕНИИ В АЗОВСКОМ МОРЕ И

ТАГАНРОГСКОМ ЗАЛИВЕ......................................96

§11. Построение сетки и диджитайзинг дна................... .96

§12. Моделирование ветровых течений........................101

ГЛАВА IV. (ЩЕНКИ СКОРОСТИ СХОДИМОСТИ СЕТОЧНЫХ УРАВНЕНИИ......114

§13. Теорема сравнения.....................................114

§14. Постановка задачи.....................................118

§15. Локально-двумерная схема..............................119

§16. Оценки скорости сходимости............................123

§17. Анализ других локально-двумерных схем.................127

ЗАКЛЮЧЕНИЕ...................................................128

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ.............................129

_ 3 -ВВЕДЕНИЕ

Актуальность темы диссертации. Большинство внутренних водоемов России и многие прибрежные районы морей являются мелководными водоемами, для которых отношение характерных размеров в горизонтальных направлениях составляет 10* и более. Но именно эти водные экосистемы испытывают наибольшую антропогенную нагрузку, связанную со строительством портов, прокладкой судоходных каналов, а также с различного рода загрязнениями. Многие из этих экосистем являются уникальными и весьма "хрупкими" (такие как Таганрогский залив и Азовское море в целом, северный Каспий и др.).

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

Среди задач водной экологии первичными являются расчеты гидрофизических параметров. Гидродинамические характеристики водоемов (поля скоростей водного потока, функции возвышения уровня) являются входными данными для расчета распространения загрязняющих веществ и развития всего живого в водной среде, начиная от фитопланктона и заканчивая млекопитающими.

Для выявления общих свойств сложных пространственных течений можно использовать численные методы [23. Но нельзя получить осмысленное решение физической задачи, связанное с нахождением полей параметров течения, если граничные условия не заданы адекватным образом. Уравнения Навье-Стокса и неразрывности представляют собой дифференциальные уравнения, описывающие удивительно широкий класс разнообразных течений, которые при изменении условий на границе могут отличаться друг от друга не только количественно, но и качественно [33.

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

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

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

Цели и задачи исследования. I. Построение гидродинамической модели, удовлетворяющей основным законам сохранения (массы, импульса, энергии) и учитывающей испарение воды и выпадение осадков, поскольку для объектов математического моделирования -Азовского моря и Таганрогского залива объем испаряющейся жидкости (преимущественно в летние месяцы) сравним с годовым речным стоком, а годовой объем выпадающих осадков уравновешивает водообмен между Азовским и Черным морями.

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

3. Дискретизация построенной модели на полученной криволинейной сетке, разработка алгоритма численного решения задачи

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

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

При выводе интегрированием по вертикальной координате трехмерных уравнений Навье-Стокса и неразрывности для вязкой несжимаемой жидкости в однородном поле тяжести во вращающейся системе отсчета получено семейство двумерных гидродинамических моделей мелководного водоема, учитывающих испарение воды и выпадение осадков С6-8], следствием уравнений баланса ишульса которых являются уравнения баланса полной механической энергии, принимающие дивергентную форму в силу уравнений баланса массы (энергетическая нейтральность) и обладающие строгой диссипатив-ностью за счет действия сил внутреннего вязкого трения, в частности, над поверхностью дна, описываемой гармонической функцией.

Многие уравнения математической физики содержат оператор Лапласа, форма которого с точностью до множителя инвариантна относительно конформных преобразований координат, минимизирующих (в случав их существования) значение функционала Дирихле в рассматриваемой области. Этот факт играет значительную роль в методах построения как регулярных, так и нерегулярных сеток. В предшествующих данному исследованию работах [9-163 была продемонстрирована высокая эффективность методов построения сеток для областей сложной формы, основанных на минимизации в расчетной области непосредственно значения функционала качества требуемой сетки, обеспечивающих независимость от выбора началь-

ного приближения. Однако, из-за чрезвычайной сложности точного решения задачи линейного поиска для функционала Дирихле в расчетной области на каждом шаге в методах 1-го порядка [9,10] не использовался функционал Дирихле. В то же время в имеющихся методах 2-го порядка [11-163 элементы матрицы Гессе [17] менялись скачкообразно, что сказывалось на работе алгоритма.

В настоящей работе за счет усовершенствований, устраняющих скачкообразный характер изменения элементов матрицы £ессе, удалось повысить эффективность имеющихся методов 2-го порядка [18,19] и, благодаря замене при точном решении задачи линейного поиска исходного функционала имеющим решение приближением, гладко передающим (до производных 2-го порядка включительно) поведение исходного функционала, удалось разработать эффективные методы 1-го порядка, направленные на получение невырожденного преобразования координат, и, используя неравенства Адамара и Коши, - расширить и обобщить классы используемых функционалов.

Использование при дискретизации по пространственным переменным конечно-элементной аппроксимации на построенной криволинейной сетке для интегралов от декартовых (а не контравариант-ных) компонент вектора скорости жидкости, а по времени - метода расщепления по физическим процессам (метода поправки к давлению) открывает дополнительную возможность - сохраняет дивергентный вид (следовательно, положительную определенность соответствующих операторов) конвективных, вязких и гравитационных слагаемых.

Для оценки решений параболических и эллиптических сеточных уравнений, главным образом, в системах координат, имеющих особенности, доказана теорема сравнения, оценивающая норму решения в равномерной метрике С выражением, содержащим норму С правой части и дискретный аналог интеграла от метрических коэффициентов монотонного сеточного оператора [20,213.

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

теплопроводности в тороидальных координатах показано, что оценки скорости сходимости при определенном расщеплении оператора (определенной локально-двумерной схеме (ЛДС)) являются неу-лучшаемыми. При доказательстве не использовалось разделение переменных и построение мажорирующих решений в явном виде (что может быть затруднительно или невозможно).

Практическая значимость. Разработанные методы построения сеток могут быть использованы для любых приложений, требующих отображения заданной области в криволинейную область более удобной формы. Построенные в результате оптимизации функционала Дирихле отображения Г-образной области (море - 50*20, залив -27*10 узлов) и прямоугольников 83*20 и 80*20 узлов на акваторию Азовского моря и Таганрогского залива максимально точно передают очертания береговой линии. По полученным на основе многолетних наблюдений 812 отметкам глубин форма дна оптимальным образом приближена интерполяционно-сглаживающей поверхностью, близкой к гармонической сеточной функции, что для песчано-илистого дна Азовского моря и Таганрогского залива является приемлемым [223. И сетки, и поверхность дна могут служить составной частью дискретных моделей других процессов в Азовском море и Таганрогском заливе. Реализованный алгоритм численного расчета ветровых течений с высокой точностью удовлетворяет закону сохранения массы. Все программы доведены до работоспособного состояния.

Достоверность полученных результатов основана на строгих математических выводах (при отмечаемых ограничениях) всех положений и максимально точном учете имеющихся входных данных, а также на качественном, а в ряде случаев и количественном соответствии результатов моделирования реально наблюдаемым течениям в Азовском море и Таганрогском заливе [233.

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

Навье-Стокса и неразрывности, дифференциального и интегрального исчисления функций многих переменных, матричной алгебры, тензорного исчисления, вариационного исчисления, методов оптимизации (нелинейного программирования), методов расщепления, метода конечных элементов, численных методов (решвдая систем линейных алгебраических уравнений, определения действительных корней полиномов), теории устойчивости разностных уравнений А.А.Самарского.

Аппробация. Результаты работы докладывались на

- XLIII Всесоюзной научной сессии, посвященной Дню радио (Москва, 1988);

- Школе-семинаре по комплексам программ математической физики (Ростов-на-Дону, 1990);

- IV Всероссийской школе молодых ученых "Численные методы механики сплошной среды" (Абрау-Дюрсо, 1992);

- IV Всесоюзном совещании "Проблемы построения сеток для решения задач математической физики" (Сунгуль (Челябинск-70), 1992);

- Областной научно-технической конференции, посвященной Дню радио (Ростов-на-Дону, 1993);

- Международной научной конференции, посвященной 100-летию со дня рождения Н.Г.Чеботарева (Казань, 1994);

- Всероссийском симпозиуме "Математическое моделирование и компьютерные технологии" (Кисловодск, 1997);

- Международной конференции "Математические модели физических процессов и их свойства" (Таганрог, 1997),

а также на научных семинарах кафедр вычислительной математики и вычислительного эксперимента, информатики, на конференциях профессорско-преподавательского состава Таганрогского государственного радиотехнического университета и на семинарах лаборатории вычислительного эксперимента вычислительного центра Ростовского государственного университета.

Автор выражает признательность своим научным руководителям

Л.С.Берштейну и А.М.Сухинову, коллективам кафедр вычислительной математики и вычислительного эксперимента и информатики Таганрогского государственного радиотехнического университета.

ГЛАВА I.

ГИДРОДИНАМИЧЕСКИЕ МОДЕЛИ МЕЛКОВОДНОГО ВОДОЕМА

В [243 содержатся теоретические основы гидромеханики, а в [25,263 - обзоры вычислительной гидродинамики и ее применений [273, направлений [283 и перспектив [293.

Гидродинамике водоемов посвящены Е30-343. Кроме этого, в следующих работах приводятся одномерная камерная модель неустановившегося течения [353, одномерные [36-393, двумерные [40423, трехмерные [433 уравнения мелкой воды и Сен-Венана [44,453, полная нелинейная система уравнений мелкой воды [463, модель многослойной мелкой воды [473 и уравнения Сен-Венана (однослойная и двухслойная модели) [483, осреднение трехмерных уравнений газовой динамики по высоте Е493, статистический вывод модели гидродинамики дисперсных сред [503, уравнения многожидкостной гидродинамики [513, модель конвекции изотермически несжимаемой жидкости [523 «двумерная Х£-модель формирования термохалинных

полей [53 У, модели диффузии вещества в движущейся среде [543, переноса токсических примесей грунтовыми водами [553, разветвленных гидравлических систем [56,573, технологии переработа полимерных материалов [583 и обрушивающихся волн [593/а в £603-

обзор моделей гидрофизики водоемов и тенденций их развития. В [613 указывается на неэквивалентность различных фор« моделей межой воды.

Различные формы уравнений Навье-Стокса можно найти в [62653, в [3,40,66,673 - способы реализации граничных условий, в [683 - аппроксимации метрических коэффициентов, а в [693 -методику анализа погрешностей и погрешности [703, обусловленные использованием криволинейных координат.

Численным методам решения уравнений Навье-Стокса посвящены [40,49,66,71-1053, в [25,26,1063 даются их обзоры, а в [1073 -обоснование метода установления. Среди них, в [40,92-943 исполь-

зуется метод конечных элементов, а в [49,66,78,87,95,98,100-105] применяется метод поправки к давлению, в т.ч. в [783 указывается на необходимость корректной постановки задач, а в [101,1023 - на необходимость сохранения производной по времени (от локального расширения потока) в уравнении Пуассона. В £53 изложен подход, дающий улучшение точности для тех задач, в которых трудно добиться выполнения уравнения неразрывности для начальных условий.

§1. Вывод уравнений гидродинамической подели

Для водоема, протяженность которого много меньше радиуса Земли, будем считать ускорение свободного падения g (с учетом центробежной силы) постоянным. В однородном поле тяжести экви-потенциаль - невозмущенная поверхность жидкости будет являться плоскостью (приближение (3-плоскости [1083). Введем прямоугольную декартову систему координат. При этом ось Oz направим противоположно направлению g из некоторой точки на невозмущенной поверхности жидкости, ось Ох - на восток, ось Оу - на север [1093. Поскольку вклад центробежной силы составляет «0,2% от вклада гравитационной силы притяжения к Земле, угол мевду вектором угловой скорости вращения Земли Q и вертикалью Oz можно считать совпадающим с широтой места -Q.

Интегрирование уравнения неразрывности

и' +v' +v>' =0

х у z

и уравнений Навье-Стокса для вязкой (в линейном приближении) несжимаемой (плотность p=const) жидкости во вращающейся с угловой скоростью

—» -»

Q=Q* (costf* J+8lnü*£), где í, J, к - единичные орты, системе координат

и;* И И / И-¡р*^' V И Н / И

V

1 +2П(и81Л$-Ж08'$), ,} -20из±ав,

Ч+ЯПисозА, )

где и=и(х,у,г,Т), у=у(х,