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

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

Плохотников Константин Эдуардович

РАЗРАБОТКА И ОБОСНОВАНИЕ НОВОГО ПОДХОДА К МЕТОДУ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ: ПРОБЛЕМА ПРЕДМЕТНОЙ МНОЖЕСТВЕННОСТИ МОДЕЛЕЙ И АНАЛИЗ АДЕКВАТНОСТИ ИХ ЦЕЛЯМ ИССЛЕДОВАНИЯ

05.13.18 — математическое моделирование, численные методы и комплексы программ

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

1 В МАЙ

005058215

Москва —2013

005058215

Работа выполнена на кафедре "Компьютерные методы физики" физического факультета Московского государственного университета имени М.В. Ломоносова

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

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

Ланеев Евгений Борисович — докт. физ.-мат. наук, профессор, Российский университет дружбы народов, профессор кафедры нелинейного анализа и оптимизации

Романовский Юрий Михайлович — докт. физ.-мат. наук, профессор, Московский государственный университет им. М.В. Ломоносова, профессор кафедры общей физики и волновых процессов физического факультета

Еленин Георгий Георгиевич — докг. физ.-мат. наук, профессор, Московский государственный университет им. М.В. Ломоносова, профессор кафедры вычислительных методов факультета ВМК

Ведущая организация: Московский государственный технический университет им. Н.Э. Баумана

Защита состоится « 14 » июня 2013 г. в 15 час. 30 мин. На заседании диссертационного совета Д 212.203.28 в Российском университете дружбы народов по адресу: Москва, ул. Орджоникидзе, дом 3.

С диссертацией можно ознакомиться в научной библиотеке Российского университета дружбы народов по адресу: 117198 Москва, ул. Миклухо-Маклая, дом 6 (отзывы на автореферат просьба направлять по указанному адресу).

Автореферат разослан «_» мая 2013 г.

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

М.Б. Фомин

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

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

Методы математического моделирования и вычислительного эксперимента как информационные технологии производства новых знаний получили в последнее время значительное развитие. Если следовать известной триаде А.А. Самарского "модель -» алгоритм -» программа"', то метод математического моделирования предполагает: изучение математических уравнений, описывающих объект исследования, постановку и проведение вычислительного эксперимента и программирование. Данные вычислительного эксперимента сравниваются с данными натурного эксперимента в том или ином смысле, при этом результаты моделирования прогнозируют результаты натурного эксперимента. Это можно отнести к прямому моделированию как к информационной технологии получения прогнозируемых, т.е. ожидаемых в результате измерений результатов. Однако можно говорить не только о прямом моделировании, в основе которого лежит идея прогноза в самом широком смысле слова, но и об обратном моделировании, как об исследовании объекта.

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

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

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

1 Самарский А.А. Проблемы использования вычислительной техники и развитие информатики// Вестн. АН СССР. 1985. №3. С.57 — 69.

2 Пытьев Ю.П. Методы математического моделирования измерительно-вычислительных систем. — М.: ФИЗМАТЛИТ, 2004. 400с.

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

В части сравнения результатов вычислительного и измерительного экспериментов исследователь формулирует критерий "приемлемого качества" прогнозирования и оценивания, например, в терминах "удовлетворительно — неудовлетворительно". В этом случае модель в результате проведенных тестов вычисления и измерения: прогнозирует измеренное значение "удовлетворительно — неудовлетворительно", оценивает модель "удовлетворительно — неудовлетворительно". При прямом вычислительном эксперименте прогнозируется измеренное значение. При обратном вычислительном эксперименте оцениваются характеристики исследуемого объекта. В той мере, в какой согласуются результаты (выводы, следствия, прогнозы, оценки) модели и измерительного эксперимента, считается, что исходная математическая модель адекватна поставленной цели исследования. Именно в связке "модель адекватна цели исследования" кроется потенциал множественной генерации моделей объекта исследования.

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

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

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

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

3 См., например: VIII Международная конференция «Идентификация систем и задачи управления» SICPRO '09. — httn://wuw.econf.mfo/SICPRQ-09

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

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

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

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

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

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

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

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

Согласно развиваемым представлениям о методе математического моделирования, результаты вычислительного эксперимента рассматриваются в качестве прогноза при прямом математическом моделировании, сравниваясь в том или ином смысле с данными натурного эксперимента. При обратном математическом моделировании производится оценивание характеристик объекта по известным экспериментальным данным. Прямое и отчасти обратное моделирование наиболее отчетливо проявились на примере построения моделей термогеометрической динамики кристалла (глава III) и турбулентности (глава IV), где, в частности, удалось добиться неплохого согласия с экспериментальными данными. На этапе как прямого, так и обратного математического моделирования разработаны различного рода численные методы, алгоритмы и программы.

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

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

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

В модели турбулентного движения (глава IV) жидкость представлялась в виде ансамбля разномасштабных дискретных жидких объектов. Для описания ансамбля выведено кинетическое уравнение, подобное уравнению Больцмана. Кинетическое

уравнение решено с помощью методов Чепмена-Энскога и Грэда. Получена замкнутая система уравнений в частных производных для описания как ламинарных, так и турбулентных течений. В части решения системы уравнений, описывающей движение турбулентной жидкости в трубе, уравнение масштабов турбулентности оказалось некорректным. Математической формой некорректности оказался член антидиффузионного типа. Некорректность была преодолена путем разделения переменных в функции масштабов с последующим интегрированием по пространству. Для численного поиска стационарных решений полученной интегро-дифференциальной системы уравнений в частных производных был разработан вычислительный алгоритм, основанный на использовании полностью неявных схем расчета с итерациями и автоматической регулировкой шага интегрирования на каждом временном слое.

Алгоритмы и комплексы программ применительно к каждой из четырех первых глав диссертации получили следующие наименования: MORPHOGENESIS*, COLLECTOR, CRYSTAL и TURBULENCE.

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

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

В главе I рассматриваемый подход к методу математического моделирования был применен к задаче описания динамики биологической формы или к проблеме морфогенеза. В данной предметной области, несмотря на огромный модельный полиморфизм, комплекс моделей обнаружить и построить не удалось. Были сформулированы следующие цели моделирования: описать динамику отдельной достаточно простой с точки зрения геометрии ткани в рамках подхода "растущий континуум" и смоделировать регенерацию ткани. Для реализации данного перечня целей была предложена и изучена новая математическая модель формообразования (морфогенеза). Уравнения модели радикально отличается от дискретных моделей на основе конеч-

* Программный комплекс MORPHOGENESIS, а также три других программных комплекса COLLECTOR, CRYSTAL и TURBULENCE зарегистрированы в Хрониках Объединенного Фонда Электронных Ресурсов "Наука и Образование" fhttr://ofernio.ru~). №02(33) февраль 2012, с. 10 — 11.

ных автоматов4 или порождающих грамматик5. В отличие от моделей, использующих уравнения в частных производных6, данная модель учитывает наличие границы растущей ткани и ставит задачу обеспечения свойства регенерации. Впервые исследовано уравнение в частных производных в тензорной форме, описывающее рост трехмерной ткани с границей. Граница при этом может быть как угодно сложно устроенной и меняться со временем. Впервые изучены уравнения, описывающие баланс некоторого вещества в пределах растущей ткани. Разработан новый алгоритм численного решения данных уравнений в одномерном случае. Построен комплекс программ MORPHOGENESIS, описывающий рост одномерной, регенерирующей ткани.

В главе II построена математическая модель коллектора электромагнитной энергии, который представляется в виде композиции множества ректенн и предназначен для преобразования различного рода шумов и электромагнитного "смога" в электрическую энергию. Данная предметная область имеет немалый модельный полиморфизм. Попытки описания ректенны, т.е. диода в совокупности с антенной, как приемника электромагнитного поля предпринимались еще с середины прошлого века7. Только в последнее время в связи с успехами микроэлектроники и нанотехноло-гий цели данного проекта приобрели отчетливое содержание8. В итоге была построена и изучена новая математическая модель коллектора. Вычислительная сложность данной задачи связана со специальным включением в цепь колебательного контура ректенны нелинейного устройства — диода и последующей композиции ректенн в решетку. Для решения уравнений модели был разработан новый алгоритм и комплекс программ COLLECTOR, который позволил провести вычислительный эксперимент и дать оценку КПД проектируемого устройства.

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

4 Нейман фонДж. Теория самовоспроизводящихся автоматов. — М.: Мир, 1971. 382с.

5 Lindenmayer A. Mathematical models for cellular interactions in development. I. Filaments with one-sided inputs// J. Theor. Biol. 1968. Vol.18, №3. P.280 — 299; Lindenmayer A. Mathematical models for cellular interactions in development. II. Simple and branching filaments with two-sided inputs// J. Theor. Biol. 1968. V.18, №3. P.300 — 315.

6 Tmng A. The chemical basis of morphogenesis// Phyl. Trans. Roy. Soc. London. 1952. Vol.237, ser.B. P.37 — 72; Wipert L. Positional information and the special pattern of cellular differentiation// J. Theor. Biol. 1969. Vol.25, №1. P.l — 47; Романовский Ю.М., Степанова H.B., Чернавский Д.С. Математическая биофизика. — М.: Наука, 1984. 304с.; Васильев В.А., Романовский Ю.М., ЯхноВ.Г. Автоволновые процессы. — М.: Наука, 1987. 300с.; Регирер С.А., Штейн А.А. Механические аспекты процессов роста, развития и перестройки биологических тканей// Итоги науки и техники. Сер. Комплексные и специальные разделы механики. Т.1. — М.: ВИНИТИ, 1985. С.З — 142.

1 СайепИ. , Welta ТА. Irreversibility and Generalized Noise// Phys. Rev. 1951. Vol.83, No.l. P.34 — 40; CallenH , BaraschM , Jackso JL Statistical Mechanism of Irreversibility// Phys. Rev. 1952. Vol.88, No.6. P. 1382 — 1386; Горелик Г.С. Некоторые применения второго закона термодинамики к электрическим флуктуа-циям// Успехи физических наук. 1952. T.44. Вып.1. С.33 — 45.

8 YmaghiHeaki , ВаШп , XhbaraNbhiko , et al. Antenna-coupled rectifying diode for IR detection// Proc. SPIE. 1996. Vol. 2882. P. 102 — 110; (Abstract: http://adsabs.harvard.edu/abs/1996SPIE.2882..102Yl.

9 Вайнштейн Б.К. Современная кристаллография. T.l. — M.: Наука, 1979. 384с.; Зоркий П.М. Симметрия молекул и кристаллических структур. — М.: Изд-во МГУ, 1986. 232с.

отношение термостата с конечным кристаллом10. В результате уравнения молекулярной динамики впервые преобразованы к форме описания ансамбля атомов образца в конфигурационном пространстве в зависимости от параметра температуры. Сформулирован алгоритм построения всей цепочки многочастичных вкладов в потенциал, представленный в виде разложения Борна-Оппенгеймера, по априорно заданной геометрии кристаллической решетки. Разработаны новые численные алгоритмы и программный комплекс CRYSTAL для расчета термической реконструкции поверхности ряда металлов. Сравнение выводов вычислительного и натурного экспериментов по термической реконструкции поверхности платины и вольфрама оказалось удовлетворительным.

В главе IV построена математическая модель турбулентности. Для описания феномена турбулентности характерен огромный модельный полиморфизм11. Множественность моделей турбулентности истолкована с точки зрения множественного толкования дискретного жидкого объекта турбулентной жидкости. Данное толкование выступает в качестве комплекса моделей. Разработана новая многомасштабная математическая модель турбулентного движения сплошной среды, т.е. модель в которой описываются дискретные жидкие объекты различных масштабов. Этапы построения многомасштабной модели включали следующие шаги. Был выведен из уравнений невязкой несжимаемой жидкости и впервые изучен бинарный потенциал, описывающий взаимодействие пары помеченных точек турбулизованной жидкости. Для данного потенциала в рамках решения задачи двух тел обнаружена особенность ветвящегося типа. На базе данного потенциала построен ансамбль разномасштабных дискретных частиц сплошной турбулентной среды. Для описания ансамбля дискретных частиц выведено кинетическое уравнение по типу уравнения Больцмана. Получена новая система многомасштабных уравнений гидродинамического типа, пригодная, как для описания ламинарных, так и турбулентных течений. В части масштабов полученная система уравнений оказалась некорректной, что потребовало использование соответствующего подхода по решению некорректно поставленных задач. Из исходной системы уравнений получена новая нелинейная многомасштабная система уравнений для описания движения турбулизованной жидкости в трубе. Разработан новый численный алгоритм решения данной системы интегро-дифференциальных уравнений в частных производных, построен комплекс программ TURBULENCE. Сравнение результатов численных расчетов движения турбулентной жидкости в трубе с экспериментальными данными оказалось удовлетворительным.

Комплексы программ MORPHOGENESIS, COLLECTOR, CRYSTAL и TURBULENCE, построенные применительно к каждой из первых четырех глав диссертации, находятся в приложенном к диссертации CD-диске в папке "Программное приложение". Все коды программ написаны в среде MATLAB.

10 Andersen В. Molecular dynamics simulations at constant pressure and/or temperature// J. Chem. Phys. 1980. Vol.72, №4. P.2384 — 2393; Nse' S. A molecular dynamics method for simulations in the canonical ensemble// Molec. Phys. 1984. Vol.52, №.2. P.255 —268.

11 Полуэмпирические теории Прандтля, Кармана, Тейлора; Колмогоров А.Н. Локальная структура турбулентности в несжимаемой жидкости при очень больших числах Рейнольдса// ДАН СССР. 1941. T.30, №4. С.299 — 303; Моисеев Н.Н. Алгоритмы развития. — М.: Наука, 1987. 304с.; МонинА.С., ЯгломА.М. Статистическая гидромеханика. 4.1. — М.: Наука, 1965. 639с.; Ктмонтович ЮЛ. Турбулентное движение и структура хаоса: Новый подход к статистической теории открытых систем. М.: Наука, 1990. 320с.; Хинце И.О. Турбулентность: ее механизм и теория. — М.: Физматгиз, 1963. 680с.

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

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

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

Результаты математического моделирования электромагнитного коллектора в главе II имеют непосредственную практическую ценность для энергетики, поскольку направлены на разработку устройства, которое, согласно замыслу, позволит преобразовать энергию, которую можно отнести к категории шума или электромагнитного "смога", в электрическую энергию. Предложенная математическая модель коллектора, а также программный комплекс COLLECTOR позволили сформулировать и подсчитать КПД устройства. Возможность априорного подсчета КПД устройства позволяет оптимизировать конструкцию коллектора в части подбора наиболее подходящей топологии ансамбля ректенн в коллекторе.

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

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

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

Апробация работы. Результаты диссертации опубликованы в 27 печатных работах [1 — 27], включая один обзор [11], две монографии [1, 13] и курс лекций [2] для студентов старших курсов. Основные положения, выводы и материалы диссертации опубликованы в монографии [1] и в курсе лекций [2].

По теме диссертации были сделаны доклады: в МГУ им. М.В. Ломоносова (семинар Ю.М. Романовского, 1978); на годовой конференции ВЦ АН СССР в г. Пущино, 1978; в Институте биологии развития им. Н.К.Кольцова (семинар А.И. Зотина, 1980); на VII Всесоюзной школе по математическому моделированию сложных биологических систем, 1980; на ряде семинаров в Гидрометцентре СССР в течение 1980— 1989 гг.; на Всесоюзной школе по теоретической физике (Калининград, 1989); на кафедре "Вычислительных методов" и в лаборатории математического моделирования в физике ф-та ВМ и К МГУ, 1989 — 2001 гг.; в Институте прикладной математики им. М.В. Келдыша РАН (семинар Г.Г. Малинецкого, 1996, 1997); на семинаре "Синергетика" (под руководством О.П. Иванова, МГУ, 2001, 2003); на семинаре "Синергетика" (под руководством Ю.Л. Климонтовича, физ. ф-т МГУ, 2001); на семинарах кафедры "Компьютерных методов физики" физического факультета МГУ, 2002; на научной конференции "Ломоносовские чтения" (сер. физики, МГУ, апрель 2003). На международной конференции "Mathematical Modeling and Computational Physics 2011", Starâ Lesnâ, High Tatra Mountains, Slovakia (July 4-8, 2011). С 2002 г. и по 2012 г. на физфаке МГУ читается спецкурс "Метод и искусство математического моделирования".

Структура и объем работы. Диссертация содержит введение, пять глав, заключение, список цитируемой литературы, а также комплексы программ, помещенные в приложенном к диссертации CD-диске. Первые четыре главы содержат введение, тематические параграфы и заключение. В последней, пятой главе диссертации в исторической ретроспективе обсуждается метод математического моделирования и феномен мультимодельности. Объем диссертации 165 страниц формата A4, 62 рисунка, 241 наименований цитируемой литературы.

Личный вклад автора в проведенное исследование. Большинство постановок проблем, предложение и исследование математических моделей, проведение соответствующих вычислительных экспериментов, а также построение всех комплексов программ MORPHOGENESIS, COLLECTOR, CRYSTAL и TURBULENCE проводилось автором самостоятельно. В процессе предложения и построения математических моделей в рамках тех или иных научных дискуссий участвовали многие исследователи из различных учреждений. Всем им мне хотелось бы выразить глубокую признательность.

СОДЕРЖАНИЕ ДИССЕРТАЦИИ

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

Кратко сформулируем основные положения нового подхода к методу математического моделирования. Более подробно они изложены во введении и в главе V диссертации.

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

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

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

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

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

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

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

ГЛАВА I. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЭЛЕМЕНТОВ МОРФОГЕНЕЗА

Морфогенез или формообразование — одно из загадочных свойств живого. В наиболее яркой форме морфогенез проявляет себя в эмбриональном развитии млекопитающих. Существует множество моделей морфогенеза, в которых толкование дискретного объекта процесса формообразования весьма разнообразно и нечто общее найти не удается. В данной предметной области к наиболее заметным именам исследователей можно отнести А.Г. Гурвича, A. Turing, Ч.М Чайлда, N. Rashensky, D'Arcy W. Thomson, R. Thorn, Дж. Фон Нейман, A. Lindenmayer, Л. Вольперта, Ю.М. Романовского, Л.В. Белоусова, С.А. Регирера, A.A. Штейна, Б.Н. Белинцева, Д.С. Чернавского, А. Gierer, Н. Meinhardt и некоторых других.

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

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

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

где / — время, Е, — фиксированная и меняющаяся области в трехмерном пространстве, л:(,) = уО,?). Предполагается существование обратного отображения г такого, что 4 = х{,) =у(1,г(Ы'))), ¿¡еЕ, хфе П(,).

Уравнение, описывающее рост ткани, имеет следующей вид:

= [(Уо + )51к + (V,(2) дЮ^ сху

где ij,k = 1,2,3 (по повторяющимся индексам в (2) предполагается суммирование); Уо — скалярный рост, Уц, Ух. — векторный рост в направлении единичного вектора и и перпендикулярно ему. Вектор и, в направлении которого осуществляется рост, можно связать с физиологическими градиентами Ч.М. Чайлда12.

Уравнение баланса некоторого вещества а в пределах растущей ткани описывается в модели выражением:

'- + -^-(Са) = ОАа + Р, С = С(Г, х) =

.......• (3)

а ох с/

где С— поле скорости движения растущей ткани, £> — коэффициент диффузии, Д — оператор Лапласа, Р — источник вещества а. В общем случае может быть несколько видов веществ, ответственных за рост ткани.

Приведем граничное условие для уравнения (3). Будем полагать, что вещество в общем случае может покидать область локализации ткани через границу области 5П(<) по закону

Сг.-Д^)

ох

= -а(А-а) 1(„, (4)

где х— единичная нормаль к поверхности 5П®, а =а((^с), А=А(^гх) — заданные на неотрицательные функции. Функция а характеризует степень проницаемости граничной поверхности для вещества а, а А — плотность вещества а вне области расположения ткани.

Система уравнений (1) — (4) при условии, что определены величины ц>, Уц, и> К Р, X> а> Л позволяет описать пространственно-временную динамику роста отдельной ткани. Численное исследование системы уравнений (1) — (4) было предпринято в простейшем одномерном случае для двух веществ: одного активатора а и одного ингибитора И роста.

После соответствующих упрощений в (1) — (4) была выведена следующая система уравнений для описания роста одномерной ткани:

12 Чайлд Ч.М. Роль организаторов в процессах развития. — М.: ИЛ, 1948. 146с.

81 1 дх ° дх2 1

дг 1 дх

а

С1(1,х) = С(1,х)-хС(1,1),

С(1,х) = у = = 'С(М);

дх

Г'Д

5А

Л

= ао(/,0), -Г'Ц,—

} &

РЬ

= Ж'.о), -'"Атг

, йс

= сгаОД);

В (5) Д,, О/, — коэффициенты диффузии активатора и ингибитора, 2Ь 5 — некоторые функции, зависящие от а и к (более подробно в диссертации), / — длина растущей ткани.

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

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

На втором этапе изучения системы уравнений (5) было предпринято численное решение нелинейной одномерной задачи. По пространству вводилась равномерная сетка с шагом 3, шаг по времени т„ выбирался автоматически на п-м временном слое, исходя из условия сходимости итераций. Рассматривалась полностью неявная разностная схема с итерациями на каждом временном слое:

(5+1) (5+1) (1+1)

(5)

-С,,

25

!* = [/ Г2 Д

(•+1) (1+1) <>+1) а.., — 2 а, +а,_,

(•) , м к1+к2[а1]1/Н1+(к6-

■*э) а,

к М(5+1)

~а1 а1 ' а„

(«1) <5 + 1) (5 + 1)

23

Д

(1 + 1) («1) А,.,-2 К 4

(5 + Ц

Д.,

(») (5 + 1) к (5) (5 + 1)

к,[а,]2 + (к6-кг) /г, —6-а, И, ;

а„

(6)

(7)

где

м

Q.. = o,

м t м ¡4 w (i) (5)

Си = £ [(0,5 о, + Y, aj + 0,5 at)S - A (i -1)<5], i = 2,..., N -1; (8)

1-2

И

С1К= 0;

(s) (s) ЛМ (s) (s)

А = [0,5 а, + £ а, + 0,5 a v ] J;

1-2

(«1)

/ _/<"> М (s+l)

-= k„(i;A-l)- I ; (9)

(s+1) (s+1) (s+1) (5+1)

a, - a <s) (I+I> (s+1) a -a («) <«+0 (»+D

^T^ = ±e I [ «1 + a2 ], -D. = la / K_]+ j.

(s+1) (s+1) (J+1) (s+1)

и - и U) (»+!) ('+!> h -h <»> <"■')

Dh = \P I [ h + h2 ], -Dh = 1 ¡3 1 [Лл._,+ hN \

о 5

(10)

Система уравнений (6) — (10) выступала в качестве вычислительного алгоритма решения исходной одномерной задачи (5) в частных производных. Уравнения (6), (7) аппроксимировали уравнения баланса активатора и ингибитора в пределах растущей ткани, уравнение (8) аппроксимировало скорость переноса, уравнение (9) — длину растущей ткани, уравнения (10) — граничные условия в точках х = 0 и х = 1 соответственно. В (6) — (10) i = 1,...Д — номера узлов сетки по пространству, s =

(0) (0) (0)

0,1,2,... — номер итераций на п-м временном слое, причем а, =а,(и),/г, = h\"\ 1 =/("> и

(s+l) (s+0 (s+l)

а<"+1> = а, , щ* = h, , /("+1) = 1 , когда в процессе счета достигалось выполнение неравенств сходимости итерационного процесса:

(s+l) (s) (S) (s+1) (J) (J) (s+1) (J) (s)

¡1 а -а\\^е1+е2\\а\\, \\ h - h |[< el + s2\\ h [|, | I - I \< e, +e21 I |, (11)

где ||*|| — некоторая норма, su s2 — неотрицательные константы, выступающие в роли абсолютной и относительной точности. Если число итераций превышало выбранное значение, шаг интегрирования по времени т„ уменьшался и расчеты начинались с начала на текущем временном слое.

Относительно систем уравнений (5) и (6), (7) в диссертации было сформулировано и доказано две теоремы. В теореме 1 доказана устойчивость в малом однородного решения системы уравнений (5) без источника активатора (ингибитора) роста. В теореме 2 доказана условная устойчивость разностных схем (6), (7). Условие выполнения теоремы 2 обеспечивалось алгоритмом автоматической регулировки шага интегрирования по времени г„.

Пример численного расчета нелинейной системы уравнений (5) в частных производных с помощью численного алгоритма (6) — (11) на базе комплекса программ MORPHOGENESIS приведен на рис.1. Согласно результатам, представленным на рис.1 можно сделать следующие выводы.

Во-первых, численно доказано наличие асимптотической стадии при / со. Это отчетливо видно на левом графике рис.1, когда длина растущей ткани выходит на некоторое предельное при t -> со значение. Во-вторых, на асимптотической стадии

распределение активатора и ингибитора по пространству в пределах растущей ткани неравномерное (графики в правой половине рис.1), что отвечает экспериментальным данным о наличии физиологических градиентов Ч.М. Чайлда. Заметим, что неоднородные по пространству распределения активатора и ингибитора являются следствием наличия у системы (5) свойства так называемой диссипативной неустойчивости.

Пространственное распределение активатора а н

Рис. 1. Динамика длины растущей ткани (левый график) и пространственное распределение активатора и ингибитора роста (правые графики)

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

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

ГЛАВА II. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КОЛЛЕКТОРА ЭЛЕКТРОМАГНИТНОЙ ЭНЕРГИИ

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

ная антенна, настроенная на ту или иную частоту. В колебательный контур, связанный с антенной включается нелинейное выпрямляющее устройство — диод. Полученный после детектирования ток, может быть использован во внешней нагрузке. Комбинация антенны и диода получила в англоязычной литературе специальное наименование — ректенна (Rectenna = Rectifier + Antenna). Если взять множество ректенн, то получим ректенную решетку, на выходе которой может получиться заметный ток.

Данный подход используется для передачи энергии на расстояние. Известно13, что в 1976 году Вильяму Брауну удалось передать СВЧ-пучком мощность 30 кВт на расстояние в 1 милю. КПД ректенны в этом эксперименте превышал 80%. Особенно амбициозными являются проекты по трансформации солнечной энергии с орбиты Земли в СВЧ излучение с последующим приемом энергии на Земле ректенной решеткой.

Наиболее интригующим выступает проект коллектора на длине волны 10,2 мкм, что соответствует максимуму излучения абсолютно черного тела при комнатной температуре 300°К. Сложность реализации данного проекта связана с разработкой диода, работающего в терагерцовом диапазоне (точнее в окрестности ~ 30 ТГц). Так, в работе8 развивается идея соединения антенны и выпрямляющего диода для детектирования электромагнитной энергии с длиной волны в окрестности инфракрасного максимума Лщи = 10,2 мкм. В качестве диода использовался барьерный диод Шотки размером 0,03 мкм, антенна выбиралась длиной 30 мкм. В работе утверждается, что на выходе связки "антенна + диод" с упомянутыми выше габаритами после воздействия С02 лазером на длине волны 10,6 мкм и мощностью 0,55 Ватт получено выпрямленное напряжение 173 нВ. Приведенные выше пространственные габариты ректенны позволяют говорить о нанотехнологической перспективе построения искомого коллектора.

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

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

1) построить приемник для каждого из источников;

2) заложить возможность съема энергии с каждого из приемников и

3) приготовить способ передачи тока в полезную нагрузку ото всех приемников коллектора.

13 Банке В.А., Лопухин В.М., Саввин В.Л. Проблемы солнечных космических электростанций// Успехи физических наук. 1977. Т. 123. Вып. 4. С.633 —655.

В качестве приемника энергии шума рассматривался обычный последовательный Л£С-колебательный контур с диодом. Съем энергии на полезную нагрузку Яы осуществляется с емкости колебательного контура.

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

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

Рис.2. Принципиальная схема коллектора

Для принципиальной схемы рис.2 были составлены уравнения, описывающие динамику токов (/^.../дг) в каждой из ректенн и в нагрузке (71оа<1). В силу того, что диод нелинейное устройство, полученная система обыкновенных дифференциальных уравнения являются нелинейной. Была разработана численная методика решения подобных систем.

Следуя законам Кирхгофа, запишем уравнения, описывающие электромагнитный коллектор. Данные уравнения представляют собой систему нелинейных обыкновенных дифференциальных уравнений N + 1 порядка для токов в коллекторе: Ц =<?!" V,/, - УгВД) -

■Лоаа — Л'У; / Ло!

(12)

где

А(7): У2У4А = 1П[1 + У5(/-А)]. (13)

Система уравнений (12), (13) записана в безразмерной форме, где уи у2, у3, у4, у5 — безразмерные коэффициенты (более подробно в тексте диссертации). Логарифм в (13) учитывает обратную вольт-амперную характеристику идеального диода. Затруднение в численном решении системы уравнений (12), (13) связано с двумя обсто-

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

Динамика мощностей шума и нагрузки

Динамика тока нагрузки

ы^Аф

Коэффициент полезного действия, СЕ

Рис.3,а. Динамика мощностей шума и нагрузки (левый Рис.3,б. Зависимость КПД коллектора (Ж=

график) и динамика тока нагрузки (правый график)

100) от величины полезной нагрузки v3

На рис.3,а приведен результат численного расчета уравнений коллектора (12), (13) состоящего из N= 100 ректенн при некоторых значениях безразмерных параметров. В левой половине рис.3,а приведены графики мгновенных значений мощностей шума Wnoise и полезной нагрузки Wload. После усреднения по отрезку времени интегрирования [0,7}] мощностей можно оценить КПД устройства по формуле: СЕ = 100%Wload/Wnoiic, где черта сверху обозначает усреднение по времени. Для данного конкретного расчета КПД коллектора оказалось на уровне 67%. На правом графике рис.3,а приведен ток в нагрузке, который является квазипостоянным. По мере роста числа ректенн N, ток нагрузки будет приближаться к некоторому постоянному значению. На рис.3,6 приведен график зависимости КПД, СЕ коллектора, который включал TV = 100 приемников шума, от величины полезной нагрузки у3. Анализ графика на рис.3,6 показывает наличие максимума КПД, что говорит о необходимости согласовывать величину полезной нагрузки с параметрами коллектора для максимально эффективного преобразования энергии шума в полезную энергию.

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

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

ГЛАВА III. МОДЕЛИРОВАНИЕ ТЕРМОГЕОМЕТРИЧЕСКОЙ ДИНАМИКИ КОНЕЧНОГО КРИСТАЛЛИЧЕСКОГО ОБРАЗЦА

Предметную область данной главы можно отнести к кристаллографической тематике в части моделирования динамики и термодинамики конечного кристаллического образца. В данной предметной области имеется заметный модельный полиморфизм. Приведем ряд имен, упомянутых в диссертации в связи с данной предметной областью: М.А. Van Hove, К. Heinz, M.K. Debe, D.A. King, К. Moller, И.П. Базаров, Э.В. Геворкян, В.В. Котенок, А. Брус, Р. Каули, J.J. Erpenbeck, W.W. Wood, W.G. Hoover, V.N. Hoover, A.J.C. Ladd, W.L. Bade, B.E. Зубкус, F. Ercolessi, E. Tossati, F. Stillinger, T. Weber, T. Shneider, J.A. Barker, B.M. Axilrod, E. Teller, T. Kihara, H. Margenau, П.М. Зоркий, P.B. Галиулин, B.C. Урусов, R. Car, H.C. Andersen, S. Nose', G. Ertl, J.M. Haile, S. Gupta, И.А. Соколова и некоторые другие.

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

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

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

Согласно методу молекулярной динамики, решаются уравнения вида: mr, =-dUI5r¡, где i=\,...,N; m, N— масса и число атомов в образце соответственно, а потенциал взят в приближении Борна-Оппенгеймера с учетом всех многочастичных вкладов:

и= Х><2)&) + '£Фт<9к) + :. + Фт(}~Ю. (14)

!-«/=;.' [-к,--

В диссертации сформулирована техника построения всей цепочки многочастичных вкладов в потенциал (14) по априорно заданной геометрии кристаллической решетки образца.

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

=> и(2) =1 (15)

У ук

где (у*:) = (у)(уАг) + + О'^О, 2У3)(г//с) = 1 •

Добавляя к С/<2) в трехчастичной форме (15) член из (14), строим трех-

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

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

В модели термогеометрической динамики конечного кристаллического образца развит двухвременной формализм, введены быстрая и медленная компоненты движения. На языке метода молекулярной динамики определены термостаты быстрой и медленной компонент. Выведены уравнения термогеометрической динамики медленной компоненты движения. В итоге были получены конструктивные уравнения для описания обратимых и необратимых фазовых переходов. Под термином "фаза" в модели понимается отдельный локальный минимум функции потенциальной энергии 11(г\,...,гц) при нулевой температуре и та конфигурация N атомов конечного образца, которая получается после повышения температуры. В модели рассмотрены бинарные фазовые переходы, т.е. переходы между парой локальных минимумов г'1' и г'2', где гм=(г1м,...,г<")),г = 1,2.

Для описания необратимого фазового перехода г(|) —> г<2) были получены следующие уравнения:

Л, =-2Г" 8 F(')=-iГX1й, (16)

тдЯ, 221. V у

где — постоянная величина имеющая размерность времени, Я = (Дь...Ду), — положение /-го атома образца в пространстве, Г— температура кристаллического образца, А2, В\ — матрица вторых и третьих производных потенциала ЩК) соответственно в локальных минимумах гт и г(2). Переход осуществлялся из в г<2) при

некоторой температуре фазового перехода 7}, которая определяется так называемой "термической силой" Т7' . Отметим, что вычисление термической силы требует заметных вычислительных мощностей, т.к. необходимо найти обратную матрицу к матрице А2, которая имеет порядок ЗNx ЗА*, где ТУ—число атомов в образце.

Для описания обратимого фазового перехода г(1)-О г<2) было получено следующее уравнение:

Д = (17)

где Л+А=], Л '

Согласно (17), обратимый переход осуществляется по схеме г1" <-» г®, где температура обратимого фазового перехода 2} = С/12/(1п йеЫ2 - 1п йеМ{), ип = и2, ии и2 — значения потенциалов в соответствующих локальных минимумах. При данном типе фазовых переходов основные затраты машинного времени идут на выбор пары подходящих локальных минимума из огромного множества возможных.

Комбинированный бинарный потенцишт Мн+сплайн

-\ «**1 0,096; /я = 4

. \ «>,09 О®)

- у Ф^Е) (сплайн) /

/

Рис.4. Комбинированный бинарный потенциал Ми + сплайн для платины

Рассмотрим применение модели к задаче описания термической реконструкции поверхностей (100) платины. Известно, что платина допускает необратимый фазовый переход на поверхности при повышении температуры от нуля до температуры фазового перехода 2} = 400°К, при этом в силу того, что решетка К является гранецентри-рованной, трехчастичный потенциал не учитывался.

Обычно при выборе бинарного потенциала ф{2) исходят из трех условий: равновесие данной решетки атомов с точки зрения однородных деформаций, соответствие модельной энергии сублимации и объемной сжимаемости образца экспериментально наблюдаемым значениям. Для учета всех трех условий выбирался потенциал в форме Ми, т.е. в форме = --¿-яг". Потенциал ф\(К) при переходил в

сплайн четвертого порядка с пятью узлами. При этом считалось, что ф (К) = 0 при Я > Яс. На рис.4 приведен искомый вид потенциала ф{1). Считалось, что К = И1' при Т — 0. Далее по мере повышения температуры на каждый атом образца действует воз-

растающая с температурой сила Fw, "подталкивающая" систему атомов не куда-то, а к переходу в г<2). Так как воздействие F{l) приводит к незначительным сдвигам, необходимо подобрать Rc в окрестности вторых соседей, т.е. Rc~4l так, чтобы при

Tf ~ 400°К система необратимо перешла в hex-фазу.

Конечный кристаллический образец платины выбирался в виде фрагмента ГЦК решетки, состоящей из трех слоев с числом атомов в каждом из них 41, 40, 41, т.е. всего 122 атома. Вся эта конструкция устанавливалась на подложку из четырех атомов так, чтобы не нарушить характерную исходную симметрию. Расчеты с помощью уравнений (16) сводились к поиску локального минимума R при различных значениях параметра температуры, т.е. строилась функция R = R(T), которая описывала движение конечного кристаллического образца в конфигурационном пространстве в зависимости от параметра температуры, когда Т е [0,7}]. Поиск локального минимума с помощью уравнений (16) осуществлялся с помощью решателей пакета MATLAB. Программные коды содержатся в приложении к диссертации комплексе программ CRYSTAL.

ф . $ . * • * - *

* . * , * .

*****

* * * *

* . * • *

* , *

* * « • «

• * " * ' * ' * '

Рис.5,а. Совместное позиционирование атомов поверхности (100)Pt до (•) и после (*) перестройки

Рис.5,б. Расположение атомов поверхности (100)W до (•) и после (*) перестройки

Результаты вычислительного эксперимента на примере моделирования реконструкции поверхности (100) платины приведены на рис.5,а. На рис.5,а приведены результаты расчета перестроенной конфигурации г(2) в части верхнего слоя (вид сверху). Расположение атомов верхнего слоя (звезды на рис.5,а) напоминает Иех-фазу перестроенной поверхности (ЮО)Рг. Точками на рис.5,а обозначены положения атомов конфигурации г(1) до перестройки.

Перейдем к изучению обратимого фазового перехода на примере моделирования термической реконструкции поверхности (100) вольфрама. Вольфрам, допускает обратимый фазовый переход на поверхности при понижении температуры ниже некоторой в окрестности 300°К, при этом использовался трехчастичный потенциал, т.к. решетка \¥ является ОЦК и треугольник ближайших соседей является равнобедренным, а не равносторонним как у решетки ГЦК.

Построим бинарный потенциал ф{2) для вольфрама. Как и в случае платины, на расстояниях в окрестности первых соседей в качестве фа) рассматривался потенциал в форме Ми. Бинарный потенциал на расстояниях в окрестности вторых соседей до-

определялся с помощью сплайна. Общий вид потенциала ф(1) для такой же, как и у 14, с тем отличием, что для число узлов сплайн-интерполяции выбиралось равным 8, а не 5. На рис.6 приведен вид бинарного потенциала ф<2) вольфрама.

Комбинированный бинарный потенциал Ми + сплайн

1 п = 8,852: т = 2,5

\ <»,№> (Ми) \ / 5 (сплайн) /

\>' / / -

V

Рис.6. Комбинированный бинарный потенциал Ми + сплайн для вольфрама

Построим трехчастичный потенциал Учитывая (15), исходим из следующего представления для ф(3):

Ф™=Т, у,Ф(2> (*/ ) + Фт Оч ,х2,х3), ;=1

(18) (19)

В (18), (19) х\, х2, хъ — длины сторон треугольника, опирающегося на произвольную тройку радиус-векторов. Пусть гь /у, гк радиус-векторы г-го, у'-го и к-то атомов кристалла, тогда хх = \r¡ - г7-|, х2 = |г, - гк\, х3 = |г,- - гк|, а V/, />„ г = 1,2,3 — некоторые параметры. В (19) функция зависящая от периметра/] треугольника, осуществляет "подрезку" трехчастичного потенциала, при этом 1, /</„;

(20)

0, />/с;

где — длина обрезания. Функция 0, в (20) строилась в виде полиномиального сплайна четвертой степени с пятью равномерно расположенными узлами интерполяции.

После ряда построений, подробности которых приводятся в тексте диссертации, параметры потенциала (18) были уточнены:

ф<3> = Мт(х1)+}/2>(х2)+^2\х3)+

(8,2184 ■ -18,6685 ■ /"' +10,5 849 ■ /3~' Ж/).

Для трехчастичного потенциала (18') проводилась проверка того, что минимальный треугольник со сторонами 1; 1; 2/л/з »1,1547 является минимумом. В частности, проверялось два условия: частные производные функции (18') по хь х2, х} должны равняться нулю, и квадратичная форма вторых производных потенциала дгф{ъ) !дхрх}, г,у =1,2,3 положительно определена. На рис.7 приведено графическое подтверждение того, что треугольник минимальных соседей ОЦК решетки является минимумов трехчастичного потенциала (18'). Поскольку трехчастичный потенциал (18') зависит от трех переменных хь хг, х3, выберем сечение х2 = 1, которое задает гиперплоскость в трехмерном пространстве х:ь х2, х3. На рис.7,а построен трехмерный графический образ потенциала (18') на гиперплоскости х2 = 1. В дополнении к рис.7,а на рис.7,б приведены линии уровня в окрестности минимального треугольника ближайших соседей, которые и доказывают наличие искомого минимума.

Потенциал

Погещиал f'*'<>. 1,'у

Рис.7,а. Трехмерное изображение потенциала (18') на гиперплоскости х2 = 1

Рис.7,б. Изображение линий уровня в окрестности минимального треугольника 1,1,2Л/3 на гиперплоскости х2 = 1

В нашу задачу входило численное моделирование термической реконструкции поверхности (100) вольфрама. Согласно ряду источников при температуре ниже 300°К, а по другим данным ниже 370°К чистая поверхность (1x1) обратимо перестраивается и образуется поверхность (V2 х л/2)Д45°, получившая название "zig-zag". Для моделирования обратимой термической перестройки поверхности W использовалось уравнение (17), где особенно важным был подбор пары локальных минимума гт, г(1) функции U, которые бы обеспечили температуру фазового перехода в окрестности экспериментального значения.

Выбирался фрагмент ОЦК решетки в виде плиты, состоящей из трех слоев с числом атомов в каждом из них 16, 9, 16, всего 41 атом. Выбранное число атомов было лимитировано имеющимися вычислительными ресурсами по поиску (путем перебора) необходимых локальных минимумов г(1) и г{2). Выбранный фрагмент ОЦК решетки помещался на подложку из пяти атомов так, чтобы не нарушать характерную симметрию атомов образца.

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

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

Г|{'+о-Г|м ас/у+") Г„ 8Г, '

где If" - U + Щ, Щ — потенциал, описывающий взаимодействие атомов образца с атомами подложки,/= 1,..., N, г/"0 = /•,(?„)> г/"+1) = r;(in+1), t„+i = tn + т„. Разностная схема в (21) является неявной, поскольку набор положений атомов г/"+1) в момент времени t„+1 не выражен в явной форме через другие известные величины. Для нахождения неизвестных положений атомов г/"4 на временном слое t„ +i определим итерационный процесс, который опишем индексом s,s = 0,1,2,... Представим набор итераций в виде:

(Sfl) (») (!)

г, =г, + &,, (22)

(0) (i+l) где г, =/;("). Подставим в (21) вместо г/"+1) -> г, в форме (22), тогда, после некоторых тождественных преобразований, получим

(23)

сг,

В дальнейшем считая, что добавка Sr, при переходе от текущей итерации к последующей мала, проведем разложение в правой части (23), тогда после некоторой перегруппировки получим следующую схему расчетного алгоритма:

. М («)

5 £/О) « , , w 8U(rt

=r{: (24)

dra,,drfik p Sr^

где Spfcui — символ Кронекера; a,p= 1,2,3; к,1= 1 Расчетную схему (24) необходимо дополнить критерием завершения итерационного процесса. Выбирая £\ и е2 достаточно малыми неотрицательными величинами, положим в качестве критерия завершения итерационного процесса (24) следующее выражение:

« «

II Зг ||<е, +г21| г ||. (25)

где ||*|| — некоторая норма. Когда верно (25) считается, что итерационный процесс

U) (»)

завершился и r,' = r, + Sr,. В случае, если итерационный процесс не сходится, уменьшается шаг интегрирования т„ и итерационный процесс повторяется.

В расчетном алгоритме (24), (25) наиболее затратной частью с точки зрения машинного времени является вычисление второй производной функции потенциальной энергии. Для ускорения процедуры решения считалось, что вторая производная потенциала U" в левой части уравнения (24) берется с нулевой итерации, т.е. при Коды для данного алгоритма написаны в среде MATLAB и содержатся в программном комплексе CRYSTAL, приложенном к диссертации в виде соответствующих файлов на CD-диске.

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

эта перестройка является обратимой. На рис.5,б приведен вид сверху на первый слой атомов вольфрама. Точками обозначены положения атомов (1х1)-поверхности (г<2)) при Т> Tß а звездочками — zig-zag-поверхности (г(1)) при Т< Tf. После построения 17 локальных минимумов г(|) и г(2), среди них была найдена такая пара, для которой температура перестройки поверхности оказалось вполне приемлемой с точки зрения экспериментальных данных и равной 7} = 315°К. Процедура подбора локальных минимумов является примером обратного моделирования, когда по данным эксперимента оценивается существование такой пары локальных минимума, которые обеспечивают нужное значение температуры фазового перехода.

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

ГЛАВА IV. МНОГОМАСШТАБНАЯ МОДЕЛЬ ТУРБУЛЕНТНОГО ДВИЖЕНИЯ ЖИДКОСТИ

Предметная область под названием турбулентность является показательной с точки зрения развиваемого подхода к методу математического моделирования. Данная предметная область обладает огромным модельным полиморфизмом. Проблема турбулентности исследовалась в трудах таких ученых, как Рейнольде, Прандгль, Карман, Тейлор, Л.Д. Ландау, А.Н. Колмогоров, A.C. Монин, A.M. Яглом, H.H. Моисеев, И.О. Хинце, Ю.Л. Климонтович, Т. Tatsumi, S.F. Edwards, В.Р. Кузнецов, Е.А. Новиков, P.M. Chung, A. Leonard, T.S. Lungren, Y.B. Pointin, Н.И. Булеев, Г.С. Глушко, А.Т. Онуфриев, В.В. Струминский, С. Кляин, И. Никурадзе, С. Улам, Ф. Харлоу, О.М. Белоцерковский, Ю.М. Давыдов, С. Херт, У.Д. Шульц, В.Е. Яницкий, В.П. Старр, Н.И. Булеев и ряда других исследователей. С точки зрения представления турбулентной жидкости в виде ансамбля дискретных жидких объектов множество моделей турбулентности можно классифицировать по тому, как толкуется в них дискретный жидкий объект. Он может быть истолкован как жидкая частица, моль, вихрь, глобула континуума и т.п. Именно эта множественность толкования дискретного жидкого объекта турбулентности положена в основу комплекса моделей турбулентности.

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

+ = (26)

dt Ал 8xldxJ8xk

где г = r(i,a), V- V(t,a) — положения в пространстве и скорости точек жидкости, а — лагранжевы координаты точек жидкости. Интеграл в (26) был получен после

решения уравнения &Р=д2(УУ])1дх1дх] относительно давления Р путем формального

обращения оператора Лапласа Д с помощью функции Грина и после подстановки полученного выражения в уравнение движения жидкости.

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

¿<«> =к<«); у(") =_ ^др^/дг™, (27)

(ЯГ

(28)

В диссертации проведено исследование потенциала взаимодействия (28) на примере изучения динамики пары фиксированных точек жидкости путем решения системы уравнений (27), которая в этом случае может быть переписана в виде: Л ="„v, =(15r(v12r12)2|r12|"7-3j'v1:!2|r12f5)r12-6/(v12r,2)]r12p5v12; r2=Vi,V2=-Vl,

где у= const >0 — некоторая константа, входящая в потенциал (28); гъ vj; г2, v2 — положение в пространстве и скорости точек 1 и 2.

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

dy \- 6х~г + (1 + 24х~3 )у2 dx ху

Переменная х в (30) характеризует расстояние между парой точек жидкости, а у — котангенс угла между относительным радиус-вектором Гц = rt - г2 и относительной скоростью V12 = vi — V2 пары точек. Качественный вид решений уравнения (30) представлен на рис.8,а.

(30)

Рис.8,а. Качественный вид решений уравнения (30)

Рис.8,6. Схема перераспределения массы, а также других аддитивных величин иллюстрирует механизм турбулентного перемешивания

Согласно рис.8,а пара точек может упасть на общий центр (движение по траектории, входящей в точку Л2), получится то, что названо в модели квазипарой. Квазипара может распасться на пару точек (движение по траектории выходящей из точки А1). Показано, что образование и распад квазипары происходит за конечное время.

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

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

Считалось, что процесс образования и распада квазипары происходит достаточно быстро. Именно в этом вся ценность понятия квазипары, т.к. она должна существовать в течение времени меньшего, чем время, затрачиваемое в среднем частицей на прохождение расстояния между двумя соударениями. Можно показать, что условие разреженности газа жидких частиц, принятое в модели, обеспечивает то, что в среднем квазипара раньше распадется, чем встретится с новой частицей. Для описания поведения ансамбля частиц вводилась функция распределения/ь где=/(f,/ï,#i), = {v\,m\,ô\,s{)\ t — время, Г\ — положение в пространстве, V] — скорость, тй[ — масса, Si — объем, Е\ — внутренняя энергия жидкой частицы. Плотность количества частиц dni в окрестности вектора имеет вид: dn\ =f\d^\, dÇ\ = dv\dm\dô\de\. В итоге оказалось возможным построить кинетическое уравнение относительно функции fif= (t,r,Ç{). Кинетическое уравнение является аналогом уравнения Больцмана:

fUv.fL-tff/ta/^ +x[(£)\i-^y\fâ(frldrdmdScjet (31)

8t or J J 4m,

где

vH^-v.J.lj'H 1;

=Ш'3 w"3+З'3). С=(£)"3[M-i«y)',3+w+ï^)"3];

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

да / ôt + tfVor = -p"'divi - v~'a,

dfH dt + wVp = -p~l (P / a)div/ + p~l(P2 / a)div/ - v~'p,

ôp/8t + div(/w) = 0,8(pw,)l St + S(_pw,wj + lf> + Pf>) / dxj = 0, (32)

d(jnT,)/ôt + div(jnT,w + qw) + ^dw, /ôxj+v'1 \nT,=0, 8(pcTJ / ôt + div(pc7> + q(m) + + /*т)Эи>,. / 8Xj -v~l^nT, = 0, где a = 1/Z(?,r), /3 = 1 ICl(t,r)\ E(t,r), Q(t,r) — средняя масса и объем жидкой частицы, а первые два уравнения системы описывают динамику масштабов турбулентности; р, w, T,, Тт — плотность, скорость, "турбулентная" и молекулярная температуры жидко-

ста. Субиндексы "/" и "т" обозначают турбулентные и молекулярные характеристики жидкости.

Условия замыкания системы уравнений (32) были получены путем приближенного решения кинетического уравнения методом Грэда. Удалось получить следующие соотношения:

I = к,а11гршт;ПЧа, J = kja-U2J32,3T,l,2V(aj3-1),

',/&,)', Р,=пТ„

qm =ik1a-1'2/S2'%,l2Va-1^(^fna',,2/32'Xl'2V(aT,), (33)

9<"» = -к;VTm, к, =ckqaV2p2'Xn-

Для замыкания системы уравнений (32) помимо (33) необходимо добавить уравнение состояния среды ЄРm(p,Tm). Величина Р,пТ , в (33) имеет смысл турбулентного давления в полной аналогии с молекулярно-кинетическим давлением идеального газа. В (32) введены молекулярные тензор Р^т) и поток энергии q<-т\ Это сделано из тех соображений, чтобы в отсутствии турбулентности, когда Т, = 0, искомые уравнения (32) совпали с системой уравнений Навье-Стокса. Такое совпадение действительно имеет место, т.к. согласно (33) при Т, = 0 — I, J, Р^, q^mf> = 0. Тем самым вопрос о соотношении молекулярных и турбулентных характеристик решен путем добавления в уравнения (32) некоторых феноменологических членов.

Система уравнений (32) допускает общий закон сохранения суммы кинетической ( 1ри>2), турбулентной (|лГ() и внутренней (рсТ„) энергий, т.е. pw1+\nT,+ pcTJ/ dt + div[(y piv2 +jnT,+ pcTJw +

(P(m) + P(") • w + W" W ] = 0.

Определим тензор of =f~') и свяжем его и поток <[mi) с хорошо известным в теории турбулентности тензором турбулентных напряжений Рейнольдса p<vJiv/j > (р= const), где <...> — операция усреднения по времени, w' = w -<и>> — пульсационная составляющая скорости и величиной p<w'T'm >, где Тт ~Тт—<Тт > — пульсационная составляющая температуры, т.е. af =р< w'j >, q(mt) = рс< w'T'm >.

В этом случае величины ¡л, и kt в (33) можно интерпретировать как коэффициенты турбулентной вязкости и теплопроводности, которые являются функциями масштабов (а, (3) и температуры турбулентности (Г,). При этом модель допускает хорошо известное в теории турбулентности приближение Буссинеска, когда тензор напряжений Рейнольдса пропорционален тензору деформаций (8 Wj/dxj) = dwjdxj + dwj/dxi - (2/3)div w-Sjj, а турбулентный поток внутренней энергии — градиенту молекулярной температуры, т.е. р< w\w'j >=-/jt{dwi IdXj) , pc<w'T'm>=-k,VTm.

Обращают на себя внимание первые два уравнения системы (31). В них был обнаружен антидиффузионный механизм, описывающий неустойчивый процесс прогрессирующего измельчения масштаба турбулентности (Р=\Ю., П — средний объем жидкой частицы), т.е. 8/3/dt + wVfi = -&\(klpU6vp)-kJ35n, где kx > 0.

Рассмотрим модельную систему уравнений (32) для случая описания турбулентного течения несжимаемой жидкости в трубе. Перейдем в цилиндрическую си-

стему координат {г, <р, £), помещая ось г в центр трубы. Положим, что давление в жидкости спадает равномерно от входа к выходу, т.е. РЛР 1 = -(АРг)/ЬР 0. где АР = Ро~Р\ — перепад давления; Р0, Р\ — давление на входе и выходе трубы соответственно, Ь — длина трубы. Согласно условию осевой симметрии и однородности характеристик по координате г, неизвестными в нашей задаче являются обратный объем жидкой частицы /?=/?(/,г), осевая компонента скорости ч> = и(/,г) и температура турбулентности Т, = Г,(Г,г). В дальнейшем удобно работать не с температурой, а с давлением турбулентности Р, =РТ,. Итак, учитывая систему уравнений (32) и условия замыкания (33), находим

81 г дг дг

д1 г дг дг рь

д1 гдг дг

(34)

г дг дг дг

где ут=цт/р — молекулярная кинематическая вязкость, а д1 =12л-(-^-)1/3 г23,39; а =_г_г1ь12/3 =011; о, =-2г=(4?У/3 з 0,18 — численные значения, найденные в ходе

построения модели турбулентного движения жидкости.

Система уравнений (34) является некорректно поставленной из-за наличия антидиффузионного члена в первом уравнении системы (34). Следуя идеям решения некорректно поставленных задач14, воспользуемся следующим способом. Представим функцию /?в разделяющихся переменных, т.е. /?(Г,г) = у{()в(г), где в(г) считается известной неотрицательной функцией, которая определена на отрезке [0,г0] и учитывает наличие ламинарного подслоя у поверхности трубы, т.е. в(г0) = 0. Подставим представление Р =ув в первое уравнение системы (34) и проинтегрируем его от 0 до г0 с весом гв{г), тогда получим следующую уже корректную систему уравнений:

2/3 , 5/3 ¿НУ 1-5. ,, „-1/3 „1/2ч

У = <*7 —Ьу , — = -— К 1 + д2Р тР'п)—] + с, 5« г дг дг

дР, 1 8 штдР

Р, -тЛ- (35)

дг г дг дг

г дг

3 г \сг )

где

1 1 11 др,.3

а = к1\Р!пв-тваг<1г1\в1г<1г, Ъ = Ч,\в%пЫг1\в\с1г, с = = сош1 >0.

рЬут

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

14 Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. — М.: Наука, 1974. 220с.

разностная схема. Итерационный процесс решения системы уравнений (35) при переходе с временного слоя п на слой п + I описывался с помощью индекса 5 = 0,1,... В итоге были получены следующие разностные схемы:

«+ч . , (а) М (I) (1+1)

= Г , (36)

(5) (5) (,Ч> (,-1)

- = ХГГ И + д й~иъ у'"3 Р(1/2) 1-г, 1'/+0,5и^'/2с7|+0,5 У (,1+0,5/ 4а

(5) «

(I) М (-1) (.-

(37)

е-и . . (5) (*) (»-и о-о

=±ГГ 0-1'3 -1/Зр1/2 Р..Ы-Р,., _

Гя 1Л+0,543^+0,5 У 1 Г,/+0,5 А2

(5) (5) (¿-1) (1-1) (5) (5) (5+1) (5 + 1)

Г-"3 ^ ■•Г""3 Ю +

(») (») (5 + 1) (5 + 1) (5) (1 + 1) (38)

^ОГ-м Г~т ^ХЛ^+^-^Г г2'3 ры +

(5) (5)

где =/(.(„),= ?„+1 = ?„ + г„,я = 0,1,... При численных расчетах

(0) (0) (0)

с помощью уравнений (36) — (38) считалось, что у =/ , = и; , Р,, =Р,, . Каждая итерация в системе уравнений (36) — (38) решалась методом прогонки. Итерации в системе уравнений (36) — (38) прекращались, когда одновременно выполнялись следующие три неравенства:

(5+1) (5) (5) (1+1) (I) (5) (5 + 1) (5) (5)

у-у%е,+е2\у |,|| м>-ы\\<е1+е2\\ч;1\\ ^-Р^^е.+е^Р,]], где ||*|| — некоторая норма. После завершения итераций при переходе со слоя п на

(5+1) (5+1) (5+1)

слой п + 1 считалось, что р = у , №*"+1> = , = .

Для решения системы уравнений (35) с помощью численного алгоритма (36) — (38) был разработан комплекс программ ТХЛШиЬЕМСЕ, который находится на СБ-диске, приложенном к диссертации. В программе шаг интегрирования по времени т„ подбирался автоматически. В процессе расчета наибольшие проблемы создавало уравнение (38). В частности необходимо было ввести проверку положительной определенности распределения значений г = 1,...Д на каждом шаге. Из-за вычислительных погрешностей значения давления турбулентности могли стать отрицательными, в этом случае значения автоматически обнулялись. В конечном счете, удалось подобрать наиболее подходящие значения параметров для решения задачи моделирования турбулентной жидкости в трубе. Коды соответствующих программ написаны в среде МАТЬАВ.

Итоги решения системы уравнений (35) с помощью численного алгоритма (36) — (38) приведены на рис.9. На рис.9,а приведен класс профилей функции 9{г), которая определяет распределение среднего обратного объема жидких частиц в трубе в зависимости от радиуса. На рис.9,б приведены три модельных профиля скорости, полученных путем решения системы уравнений (35) с помощью численного алгоритма (36) — (38) применительно к задаче описания движения турбулентной жидкости в трубе при различных значениях перепада давления (с). Каждый профиль скорости

сравнивался с экспериментальными значениями, отмеченными на рис.9,б звездочками. Рис.9,б демонстрирует неплохое соответствие модельных и экспериментальных распределений скорости в зависимости от радиуса трубы. Заметные отличия модельных и экспериментальных значений скорости наблюдаются в пристеночной области.

Рис.9,а. Типичный вид графика функции 0(г)

Рис.9,б. Сравнение модельных профилей скорости в трубе с экспериментальными значениями13

Критерием адекватности многомасштабной модели турбулентного движения жидкости цели исследования выступили: итоговая система уравнений в частных производных (32), (33), а также пример сравнения данных вычислительного и натурного экспериментов течения турбулентной жидкости в трубе (рис.9,б). Сравнение данных модели вычислительного эксперимента на примере движения жидкости в трубе с экспериментальными данными оказалось удовлетворительным. Модель турбулентного движения жидкости удовлетворительно верифицирована на примере построения и численного решения уравнений, описывающих движение турбулизованной жидкости в трубе. В части реконструкции профиля распределения обратного объема жидких дискретных частиц (рис.9,а) решалась обратная задача моделирования, т.е. по экспериментальному значению профиля скорости был восстановлен профиль распределения масштабов турбулентности.

ГЛАВА V. МЕТОД МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ И ПРОБЛЕМА МУЛЬТИМО ДЕЛЬНОСТИ

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

Суммирование современной научной методологии и системного подхода, позволяет в любой модели фрагмента предметной области выделить следующие три характерных признака: 1) определяется фрагмент предметной области как исходная целостность, 2) выделяется дискретный элемент фрагмента и 3) осуществляется синтез

15 Никурадзе И. Закономерности турбулентного движения жидкостей в гладких трубах// Проблемы турбулентности. — М.-Л.: ОНТИ, 1936. С.75 — 150.

ансамбля дискретных элементов фрагмента в новую целостность. На языке системного подхода, рядом авторов16 модель фрагмента представляется в виде схематической цепочки: целое -» ансамбль элементов -> целое'. Математической модель становится в той мере, в какой используется язык математики для описания элемента части предметной области, его целостности, а также переходов "целое -> ансамбль элементов" и "ансамбль элементов —> целое"'. Появление множества, вообще говоря, различных математических моделей, описывающих разные части данной предметной области, зачастую, можно связать с тем, что переход "целое —> ансамбль элементов —> целое"' процедура весьма неоднозначная.

Изучение истории развития метода математического моделирования на примере различных предметных областей показало, что методология математического моделирования эволюционирует в своем развитии по некоторому стандартному трех-этапному образцу. В начале, до всякого анализа, предметная область выступает как неразделенное целое, цели моделирования не определены. На этапе I имеет место од-на-единственная математическая модель предметной области, которую в силу ее единственности зачастую догматизируют (пример: классическая механика до создания специальной теории относительности). Для первого этапа характерно отчетливое разделение понятий "модель", "натурный эксперимент" и сравнение данных модели с данными натурного эксперимента, т.е. на этом этапе вопрос об адекватности модели объекту исследования не вызывает затруднений. На этапе II у первоначально доминирующей модели появляются конкуренты. Предметная область расщепляется на несколько фрагментов, для описания которых строятся, вообще говоря, разные математические модели. Для второго этапа взаимоотношения между моделями и так называемым натурным экспериментом значительно усложняются и опосредуются. Вопрос об адекватности моделей фрагментам предметной области трансформируется в вопрос об адекватности цели исследования, которую формулирует исследователь. Второй этап переходит в третий этап с появлением комплекса моделей, который сдерживает процесс безудержного производства все новых и новых моделей фрагментов предметной области, отличающихся лишь по форме. Наличие комплекса моделей позволяет охарактеризовать имеющиеся модели с единых позиций. Этап зрелости моделирования той или иной предметной области характеризуется этапом III, когда комплекс моделей данной предметной области считается созданным. Для третьего этапа характерно наблюдение и моделирование предметной области во всей своей полноте.

Определим, что понимается под математическим моделированием предметной области и под математической моделью части предметной области.

ОПРЕДЕЛЕНИЕ №1. Под математическим моделированием некоторой предметной области понимается специфически человеческая деятельность, развертывающаяся поэтапно (этап I, этап II, этап III) и приводящая к созданию комплекса моделей, который является системной сборкой отдельных моделей фрагментов предметной области.

16 Морозов К.Е. Математическое моделирование в научном познании. — М.: Мысль, 1969. 212с.; Новик КБ. О моделировании сложных систем. — M.i Мысль, 1965. 335с.j Штпофф В.А. Роль моделей в позне-нии. — Л.: Изд-во ЛГУ, 1963.128с.

ОПРЕДЕЛЕНИЕ №2. Математической моделью части предметной области является такая конструкция, которая содержит:

1) специфические для данного фрагмента предметной области определения: части предметной области в целом, дискретного элемента фрагмента предметной области и, соответственно, ансамбля элементов;

2) математический аппарат, описывающий: отдельный элемент фрагмента предметной области, ансамбль элементов, а также переходы между ними.

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

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

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

Сформулируем критерий адекватности модели цели исследования в части прямого и обратного моделирования. Критерий адекватности модели включает постановку и проведение вычислительного и так называемого натурного экспериментов. Натурный эксперимент в современном смысле слова это, как правило, сложная комбинация объекта измерения, среды, а также измерительной и вычислительной компонент в рамках единого измерительно-вычислительного комплекса*. Можно говорить о модели измерительного эксперимента и о ее качестве в части адекватности цели исследования, обусловившей проведение измерительного эксперимента. Для подтверждения качества адекватности модели цели исследования, как правило, проводится тестовый вычислительный эксперимент и тестовый измерительный эксперимент. В режиме тестирования модель должна предоставить: возможность компьютерного моделирования прогнозирования и оценивания (интерпретации измерения), причем с максимальной гарантированной точностью, а также обеспечить эмпирическую проверку качества прогнозирования и оценивания.

В рамках используемой терминологии исследователь может сформулировать критерий "приемлемого качества" прогнозирования и оценивания, например, в терминах "удовлетворительно — неудовлетворительно". В этом случае модель в результате проведенных тестов вычисления и измерения: прогнозирует измеренное значе-

* Здесь и далее, включая определение адекватности, учтены замечания Ю.П. Пытьева

ние "удовлетворительно — неудовлетворительно", оценивает модель "удовлетворительно — неудовлетворительно". При прямом вычислительном эксперименте прогнозируется измеренное значение. При обратном вычислительном эксперименте оцениваются характеристики исследуемого объекта. В той мере, в какой согласуются результаты (выводы, следствия, прогнозы, оценки) модели измерительного эксперимента, считается, что и исходная математическая модель адекватна цели исследования.

Критерий адекватности. Математическая модель фрагмента предметной области:

1) адекватна цели исследования, когда и прогнозирование н оценивание удовлетворительны;

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

3) адекватна цели исследования в части оценивания, когда модель интерпретации измерений удовлетворительна, а модель прогнозирования измерений неудовлетворительна;

4) не адекватна, как в части прогнозирования, так и в части оценивания.

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

1) добиться удовлетворительного прогнозирования (прямое моделирование);

2) обеспечить удовлетворительное оценивание (обратное моделирование);

3) разработать комплекс моделей (метамоделирование).

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

Таблица Л контексте оп »1. Классификация частных математических моделей диссертации в ределения Д°1 математического моделирования предметной области

1 2 3

Глава ДУЙ Этап ДоДо Комплекс моделей

Глава I Этап II отсутствует

Глава II Этап II отсутствует

Глава III Этап II отсутствует

Глава IV Этап III Множественное толкование дискретного объекта турбулентности

В таблице №2 приведена сводная классификация частных математических моделей диссертации в контексте определения №2 математической модели фрагмента предметной области. Согласно определению математической модели фрагмента предметной области, представленному выше, в таблице №2 во втором столбце приводится характеристика исходной целостности, в третьем столбце — дискретного эле-

мента фрагмента предметной области, а в четвертом столбце — новая характеристика предметной области после построения модели.

Таблица №2. Классификация частных математических моделей диссертации в контексте определения №2 математической модели фрагмента предметной области

1 2 3 4 5

Глава Предметная область Дискретный элемент фрагмента предметной области Предметная область после построения математической модели Комплекс моделей

Глава I Морфогенез, формообразование Элементарная ткань в форме параллелепипеда, цилиндра Регенерирующая форма, как растущий континуум отсутствует

Глава II Коллектор некогерентной распределенной в пространстве электромагнитной энергии Ректенна Система нелинейных дифференциальных уравнений, описывающих работу коллектора отсутствует

Глава III Проблема дальнего порядка в кристаллографии Отдельные атомы Дальний порядок — следствие специального вида потенциала взаимодействия в форме суммы многочастичных вкладов отсутствует

Глава IV Проблема турбулентности Дискретный элемент сплошной среды Единая система уравнений в частных производных для описания, как ламинарных, так и турбулентных течений Множественное толкование дискретного объекта турбулентности

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

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

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

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

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

ЗАКЛЮЧЕНИЕ

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

РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

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

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

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

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

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

5. В предметной области, посвященной турбулентности, определен комплекс моделей, сформулированный с точки зрения множественной интерпретации дискретного жидкого объекта турбулентной жидкости. На базе комплекса моделей турбулентности предложена и изучена многомасштабная модель турбулентности, в которой выведен и впервые изучен бинарный потенциал, описывающий взаимодействие пары фиксированных точек турбу-лизованной сплошной среды. Обнаружена особенность ветвящегося типа данного потенциала. Ансамбль дискретных объектов описан с помощью уравнения Больцмана, из которого получены уравнения гидродинамического типа, описывающие многомасштабное движение турбулентной жидкости. Выведены нелинейные многомасштабные уравнения для описания движения турбулентной жидкости в трубе. Для решения уравнений, описывающих движение турбулентной жидкости в трубе, разработан вычислительный алгоритм и комплекс программ TURBULENCE, позволивший осуществить сравнение данных вычислительного и натурного экспериментов. Сравнение данных модели вычислительного эксперимента на примере движения жидкости в трубе с экспериментальными данными оказалось удовлетворительным. Параметры модели турбулентного движения жидкости удовлетворительно верифицированы на примере построения и численного решения уравнений, описывающих движение турбулизованной жидкости в трубе.

РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ К.Э. ПЛОХОТНИКОВА ОПУБЛИКОВАНЫ В СЛЕДУЮЩИХ ИЗДАНИЯХ*

1. Плохотников К.Э. (монография) Математическое моделирование и вычислительный эксперимент. Методология и практика. — М.: Эдиториал УРСС, 2003. 282с.

2. Плохотников К.Э. Метод и искусство математического моделирования: курс лекций. — М.: Флинта, 2012. 518с. — ISBN 978-5-9765-1541-3.

3. Плохотников К.Э. Математическое моделирование и вычислительный эксперимент. Методология и практика. Научная конференция "Ломоносовские чтения". Секция физики. 18-25 апреля 2003 г. Сборник расширенных тезисов докладов. — М.: Физический факультет МГУ, 2003. С.27 —29.

4. Плохотников К.Э. Математическое моделирование и вычислительный эксперимент: методология и практика// Интеллектуальные системы. 2009. Т.13. Вып.1-4. С.5 — 32.

5. Плохотников К.Э. Модель роста одномерной ткани// Биофизика. 1982. Т.27. Вып.4. С.689 — 693.

* Полужирным шрифтом выделены публикации автора диссертации в рецензируемых изданиях из списка ВАК

6. Плохотников К.Э. Морфогенез клеточных пластов/ Математическая биология развития. — М.: Наука, 1982. С. 155 —159.

7. Плохотников К.Э. Коллектор некогерентной распределенной в пространстве электромагнитной энергии// Математическое моделирование. 2009. Т.21, №12. С.35 — 46.

8. Plbhtniko К Collector of Noncoherent Electromagnetic Energy Distributed in Space// Mathematical Models and Computer Simulations. 2010. Vol.2, No.4. P.504 — 513.

9. Плохотников К.Э. Как возможен конечный кристалл при нулевой температуре (энергети-ко-геометрическнй аспект)?// ДАН СССР. 1991. Т.320, Л»4. С.877 — 881.

10. Плохотников К.Э. Термогеометрическая динамика конечного кристаллического образца// ДАН. 1992. Т.326, Л»6. С.994 — 998.

11. Плохотников К.Э. (обзор) Термогеометрическая динамика конечного кристаллического образца// Математическое моделирование. 1993. Т.5, ЛИ. С.З — 31.

12. Plokhotnikov К.Е. Mathematical Modeling of Thermal Restructuring of the Platinum Surface// Вестннк Российского университета дружбы народов. Серия: Математика, информатика, физика. 2011, Л»3, с.69 — 75.

13. Плохотников К.Э. (монография) Математическое моделирование. Экзистенциальный аспект. — М.: Изд-воМГУ, 1993. 224с.

14. Плохотников К.Э. К вопросу о моделировании турбулентного движения. — Препринт ИПМ АН СССР, №7. — М., 1980. 26с.

15. Плохотников К.Э. Модель турбулентности, изотропной и однородной по моменту импульса квазипары. — Препринт ИПМ АН СССР, №8. — М., 1980. 29с.

16. Плохотников К.Э. Многомасштабная модель турбулентности. — Препринт ИПМ АН СССР, №64, —М., 1980. 27с.

17. Плохотников К.Э. Об одной математической модели турбулентного движения жидкости// Тр. Гидрометцентра СССР. Вып.248. —JI.: Гидрометеоиздат, 1982. С.52 — 66.

18. Плохотников К.Э. Об одной математической модели турбулентного движения жидкости// Тр. Гидрометцентра СССР. Вып.273. —JL: Гидрометеоиздат, 1986. С.15 —48.

19. Плохотников К.Э. Об одной математической модели турбулентного движения жидкости// Математическое моделирование. Современные проблемы математической физики и вычислительной математики. — М.: Наука, 1989. С.270 — 284.

20. Плохотников К.Э. Об одной математической модели турбулентного движения жидкости// ДАН СССР. 1988. Т.301, „\»4. С.805 — 809.

21. Plokhotnikov К.Е. Numerical Description of the Fluid Pipe Motion with Multiscale Turbulence Model// Вестпик Российского университета дружбы народов. Серия: Математика, информатика, физика. 2011. JV»4. С.107 — 112.

22. Плохотников К.Э. Общая циркуляция атмосферы: синтез дискретного подхода и лагранжевого способа описания сплошной среды// Изв. АН СССР. Сер. Физика атмосферы и океана. 1985. Т.21, №10. С.1110 —1111.

23. Плохотников К.Э. Математическая модель, объединяющая дискретные представления и лагранжев способ в описании сплошной среды// ДАН СССР. 1987. Т.294, -V» 1. С.75 — 79.

24. Плохотников К.Э. Общая циркуляция атмосферы: синтез дискретного подхода и лагранжевого способа описания сплошной среды. Деп. в ВИНИТИ 3.06.1985, №3827-85. 31с.

25. Плохотников К.Э. Исследование движения атмосферы в терминах ансамбля воздушных масс// Тр. Гидрометцентра СССР. Вып.273. —Л.: Гидрометеоиздат, 1986. С.49 — 79.

26. Плохотников К.Э. Общая циркуляция атмосферы: синтез дискретного подхода и лагранжевого способа описания сплошной среды// Тр. Гидрометцентра СССР. Вып.286. — Л.: Гидрометеоиздат, 1987. С.64 —84.

27. Плохотников К.Э. Программный комплекс MORPHOGENESIS для моделирования элементов морфогенеза; Программный комплекс COLLECTOR для математического моделирования коллектора электромагнитной энергии; Программный комплекс CRYSTAL, позволяющий моделировать термическую реконструкцию поверхности ряда металлов; Программный комплекс TURBULENCE для моделирования движения турбулентной жидкости в трубе с помощью многомасштабной модели турбулентности// Хроники Объединенного Фонда Электронных Ресурсов "Наука и Образование", №02(33) февраль 2012, с.10 — 11.

Краткая аннотация диссертации К.Э. Плохотникова "Разработка и обоснование нового подхода к методу математического моделирования: проблема предметной множественности моделей и анализ адекватности их целям исследования", выдвинутой на соискание ученой степени доктора физико-математических наук по специальности 05.13.18 — математическое моделирование, численные методы и комплексы программ.

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

Brief summary of the dissertation by K.E. Plokhotnikov "Development and justification of a new approach to the method of mathematical modeling: the problem of subject plurality of models and analysis of the adequacy of their research goals" putted forward for an academic degree of doctor of physical-mathematical sciences obtaining, specialty 05.13.18 — mathematical modeling, numerical methods and software.

Dissertation is devoted to the development of a new approach to the method of mathematical modeling in part of problem of subject plurality of models and analysis of the adequacy of their research goals. In a new approach to the method of mathematical modeling a multi-model is considered in terms of the adequacy of models variety of general and specific objectives of the study. This approach is demonstrated by the examples of constructing of four mathematical models in four subject areas: morphogenesis, the electromagnetic collector, thermal geometry dynamics of the final crystal sample and turbulence.

Напечатано с готового оригинал-макета

Подписано в печать 18.03.2013 г. Формат 60x90 1/16. Усл.печ.л. 2,0. Тираж 100 экз. Заказ 2698.

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к. Тел. 8(495)939-3890/91. Тел./факс 8(495)939-3891.

Отпечатано в ППП «Типография «Наука» 121099, Москва, Шубинский пер., 6

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

имени М.В.ЛОМОНОСОВА

ФИЗИЧЕСКИЙ ФАКУЛЬТЕТ

0520іо5іі)ІЗ На правах рукописи

Плохотников Константин Эдуардович

РАЗРАБОТКА И ОБОСНОВАНИЕ НОВОГО ПОДХОДА К МЕТОДУ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ: ПРОБЛЕМА ПРЕДМЕТНОЙ МНОЖЕСТВЕННОСТИ МОДЕЛЕЙ И АНАЛИЗ АДЕКВАТНОСТИ ИХ ЦЕЛЯМ ИССЛЕДОВАНИЯ

05.13.18 — математическое моделирование, численные методы и комплексы программ

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

Москва —2013

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ..................................................................................................4

ГЛАВА 1......................................................................................................15

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

МОРФОГЕНЕЗА...............................................................................................15

§ 1. Введение............................................................................................................15

§2. Рост отдельной ткани.......................................................................................16

§3. Баланс вещества в пределах растущей ткани................................................19

§4. Одномерное приближение...............................................................................20

§5. Рост одномерной ткани. Вычислительный эксперимент.............................27

§6. Моделирование роста трех связанных одномерных тканей........................36

§7. Заключение.......................................................................................................41

ГЛАВА II.....................................................................................................43

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КОЛЛЕКТОРА ЭЛЕКТРОМАГНИТНОЙ ЭНЕРГИИ..............................................................43

§1. Введение............................................................................................................43

§2. Постановка задачи............................................................................................45

§3. Приемник шума................................................................................................46

§4. Численное решение уравнений приемника шума.........................................48

§5. Коллектор электромагнитной энергии...........................................................53

§6. Численное решение уравнений коллектора...................................................54

§7. Источники энергии, отличающиеся от белого шума....................................58

§8. Заключение.......................................................................................................61

ГЛАВА III...................................................................................................63

МОДЕЛИРОВАНИЕ ТЕРМОГЕОМЕТРИЧЕСКОЙ ДИНАМИКИ КОНЕЧНОГО КРИСТАЛЛИЧЕСКОГО ОБРАЗЦА.....................................63

§1. Введение............................................................................................................63

§2. Как возможен конечный кристалл при нулевой температуре?...................65

§3. Двухвременной формализм.............................................................................69

§4. Окрестность нулевой температуры................................................................76

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

реконструкции поверхности (100)Р1................................................................................80

§6. Вычислительный эксперимент на примере моделирования

реконструкции поверхности (100)"№................................................................................89

§7. Заключение.....................................................................................................103

ГЛАВА IV.................................................................................................105

МНОГОМАСШТАБНАЯ МОДЕЛЬ ТУРБУЛЕНТНОГО

ДВИЖЕНИЯ ЖИДКОСТИ............................................................................105

§1. Введение..........................................................................................................105

§2. Исследование потенциала взаимодействия.................................................109

§3. Вывод и решение основного кинетического уравнения.............................115

§4. Исследование вопроса об измеряемости......................................................121

§5. Пример расчета турбулентного течения жидкости в трубе.......................128

§6. Заключение.....................................................................................................136

ГЛАВА V...................................................................................................138

МЕТОД МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ И ПРОБЛЕМА МУЛЬТИМО ДЕЛЬНОСТИ....................................................138

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

РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ........................................................153

РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ К.Э. ПЛОХОТНИКОВА ОПУБЛИКОВАНЫ В СЛЕДУЮЩИХ ИЗДАНИЯХ..................................154

ЛИТЕРАТУРА.........................................................................................156

ВВЕДЕНИЕ

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

Методы математического моделирования и вычислительного эксперимента как информационные технологии производства новых знаний получили в последнее время значительное развитие. Если следовать известной триаде A.A. Самарского "модель —» алгоритм -» программа" [1], то метод математического моделирования предполагает: изучение математических уравнений, описывающих объект исследования, постановку и проведение вычислительного эксперимента и программирование. Данные вычислительного эксперимента сравниваются с данными натурного эксперимента в том или ином смысле, при этом результаты моделирования прогнозируют результаты натурного эксперимента. Это можно отнести к прямому моделированию как к информационной технологии получения прогнозируемых, т.е. ожидаемых в результате измерений результатов. Однако можно говорить не только о прямом моделировании, в основе которого лежит идея прогноза в самом широком смысле слова, но и об обратном моделировании, как об исследовании объекта.

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

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

Обратное моделирование, которое нацелено на изучение внутренних характеристик объекта, имеет непосредственное отношение к проблеме множественности, мультимодельности математических моделей, используемых для описания объекта. И для прямого, и для обратного математического моделирования на этапе верификации модели готовится один или несколько тестовых вычислительных и измерительных экспериментов. Полученные результаты сравниваются в рамках выбранной или приемлемой точности, что дает возможность оценить или скорректировать структуру и параметры модели. Когда этап верификации модели завершен, определяются условия для проведения нового эксперимента и по результатам моделирования прогнозируется, какой должен быть результат. Для обратного моделирования на этапе отбора моделей из потенциального множества характерны процедуры структурной и параметрической идентификации*.

В части сравнения результатов вычислительного и измерительного экспериментов исследователь формулирует критерий "приемлемого качества" прогнозирования и оценивания, например, в терминах "удовлетворительно — неудовлетворительно". В этом случае модель в результате проведенных тестов вычисления и измерения: прогнозирует измеренное значение "удовлетворительно — неудовлетворительно", оценивает модель "удовлетворительно — неудовлетворительно". При прямом вычислительном эксперименте прогнозируется измеренное значение. При обратном вычислительном эксперименте оцениваются характеристики исследуемого объекта. В той мере, в какой согласуются результаты (выводы, следствия, прогнозы, оценки) модели и измерительного эксперимента, считается, что исходная математическая модель адекватна поставленной цели исследования. Именно в связке "модель адекватна цели исследования" кроется потенциал множественной генерации моделей объекта исследования.

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

В диссертации возможности рассматриваемого подхода иллюстрируются на ряде примеров математического моделирования: элементов морфогенеза, электромагнитного коллектора, конечного кристалла и турбулентно-

* См., например: VIII Международная конференция «Идентификация систем и задачи управления» SICPRO '09. — http://www.econf.info/SlCPRQ-09

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

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

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

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

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

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

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

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

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

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

Методы исследования. Декларируется общая для всех предложенных моделей цель: уравнения математической модели должны допускать преобразование в конструктивную форму, допускающую вычислительный эксперимент. В модели морфогенеза (глава I) уравнения становятся вполне конструктивными после перехода к одномерному случаю. В модели электромагнитного коллектора (глава II) уравнения допускают вычислительный эксперимент после специального включения диода в электрическую цепь. В модели термогеометрической динамики конечного кристаллического образца (глава III) пройден длинный путь теоретических выкладок прежде, чем удалось преобразовать уравнения к форме, допускающей численное описание обратимых и необратимых фазовых переходов. Наконец, в многомасштабной модели турбулентности (глава IV) от статистического описания ансамбля дискретных жидких объектов к уравнениям в частных производных гидродинамического типа и далее уже к конструктивным уравнениям, описывающим движение турбулентной жидкости в трубе.

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

объекта по известным экспериментальным данным. Прямое и отчасти обратное моделирование наиболее отчетливо проявились на примере построения моделей термогеометрической динамики кристалла (глава III) и турбулентности (глава IV), где, в частности, удалось добиться неплохого согласия с экспериментальными данными. На этапе как прямого, так и обратного математического моделирования разработаны различного рода численные методы, алгоритмы и программы.

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

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

В модели термогеометрической динамики конечного кристаллического образца (глава III) вычислительная задача свелась к поиску локального минимума функции потенциал