«МОДЕЛИРОВАНИЕ ПОДВОДНОГО ВЗРЫВА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ ДИССЕРТАЦИЯ на соискание учёной степени кандидата физико-математических наук ...»

Пятая глава рассматривает изменение формы ударной волны со временем.

Данное явление анализируется при помощи подхода, применяемого для солитонных явлений. Приводится вывод уравнения для двумерных ударных волн, обобщающего уравнение КдВ для одномерной среды. В полярных координатах оно имеет вид + + (, ) = 0.

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

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

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

Основой компьютерного моделирования является математическая модель исследуемого физического явления, записанная в виде дифференциальных или интегро-дифференциальных уравнений; при этом также формулируются допущения и границы применимости модели. Затем определяющие уравнения подвергаются дискретной алгебраической аппроксимации, для того чтобы можно было получить численное решение при помощи компьютера. Метод дискретизации непосредственно зависит от каждой конкретной задачи и может включать в себя методы конечных разностей, методы конечных и граничных элементов, а также методы частиц (разновидностью которых и является метод молекулярной динамики). В конечном итоге производится вычислительный эксперимент путем реализации в компьютерной программе, на выбранном языке программирования, определенных численных методов и алгоритмов решения полученной системы алгебраических уравнений. Последним шагом является непосредственная обработка полученных результатов и их визуализация. Схематично данная последовательность действий представлена на Рис. 1.1.

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

Рисунок 1.1: Схема вычислительного эксперимента 1.2 Общие сведения о подводном взрыве Взрывом называется любой процесс быстрого превращения одного вида энергии в другой. Данный процесс характеризуется высокоскоростными изменениями в физическом и химическом состоянии среды, сопровождается практически мгновенным выделением тепла, увеличением давления и температуры, возникновением больших напряжений в среде, разрушением материала и неустановившимся движением среды. Химические среды, способные к превращениям такого рода, называют взрывчатыми веществами. Процессы, происходящие со скоростями порядка 1 км/c и выше, называют взрывными. Ввиду сложности механизма энергетического превращения и его различности в каждом конкретном случае, при постановке задачи о действии данного процесса часто используются различные приближения и модельные представления [52], поскольку в общем случае распределение физических величин на поверхностях раздела произвольно. Именно существование поверхности раздела двух сред оказывается несовместимым с законами сохранения массы, количества движения и энергии. Происходит мгновенный распад начальной поверхности с образованием новых поверхностей. На каждой из этих поверхностей разрыва такие гидродинамические переменные, как давление, плотность, температура и скорости частиц, меняются скачкообразно. Нестационарная поверхность сильного разрыва, которая распространяется по окружающей среде, носит название «фронт ударной волны»;

поверхность, разделяющая продукты взрыва и среду, также является поверхностью сильного разрыва (ее называют поверхностью «газового пузыря»). Таким образом, теоретические исследования касались, в основном, изучения неустановившегося движения жидкости между граничными поверхностями: фронтом ударной волны, поверхностью газового пузыря и продуктами детонации взрывчатых веществ [53–55]. В данных исследованиях иногда выделяют такие подзадачи как образование ударной волны [56] и развитие образовавшегося при этом пузыря [57].

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

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

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

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

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

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

1.3 Краткое описание предыдущих исследований Задача о подводном взрыве являлась темой многочисленных исследований.

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

В исследовании подводных взрывов, согласно [3], принято выделять такие проблемы как ударные волны и динамика полости с продуктами детонации. При этом предполагается, что рассматриваемая среда безгранична [58]. Для решения этих проблем принято определять распространение ударной волны при помощи какого-либо приближенного метода и упрощать последующий анализ динамики пузыря, рассматривая воду как несжимаемую жидкость, а пузырь как однородный газовый шар [1]. Данная модель идеальной несжимаемой жидкости является наиболее распространенным и простым приближением многих задач гидродинамики взрыва и, как известно, предполагает постоянство плотности жидкости = и позволяет в рамках законов сохранения исследовать ряд важных характеристик взрывного процесса. Однако при значительных давлениях, возникающих при взрывах, жидкости являются довольно сильно сжимаемыми. При давлениях в 10 000 атм. вода сжимается более чем на 20 процентов. Кроме того, предположение о несжимаемости среды не дает возможности учитывать рассеивание энергии и тем самым исключает из анализа диссипативные процессы.

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

Что касается натурных экспериментов, то в литературе встречаются достаточно скудные сведения о них. В качестве примера, можно привести эксперимент мощного подводного взрыва, который был произведен 14 мая 1955 года и состоял в детонации 30-килотонного ядерного устройства (“Вигвам”), расположенного на глубине 610 м в водах Тихого океана. Целью операции “Вигвам” было определение уязвимости подводных лодок к глубоководному взрыву, и возможность использования таких взрывов в реальной боевой ситуации. В ходе данного эксперимента проводились подводные измерения давления в зависимости от времени в различных местах, регистрировались сигналы ударной волны, волны давления, фиксировались визуальные эффекты на поверхности воды. После этого была проведена значительная экспериментальная работа, как при помощи лабораторных масштабных исследований, так и при помощи сравнительно небольших полевых испытаний с взрывчатыми веществами [59]. Эта работа привела к развитию полуэмпирических моделей для описания, например, поведения газового пузыря при подводном взрыве. Однако данные таких моделей не дают полной информации относительно действительной картины подводного взрыва.

1.4 Возникновение волн при подводном взрыве Тема возникновения волн, в отличие от других явлений, происходящих при подводном взрыве (ударная волна, движение газового пузыря, пульсации давления, влияние свободной поверхности и другие процессы, сопровождающие взрыв) не так популярна. Тем не менее, данную тему нельзя игнорировать, поскольку взрывы около поверхности могут порождать волны значительной величины. Теоретический подход к данному вопросу, главным образом, основан на гидродинамической феноменологии зарождения и распространения волн [41]. В отличие от такого подхода в данной работе при компьютерном моделировании был использован метод молекулярной динамики, который позволяет получить результаты без осреднения микро-характеристик. Можно отметить, что данный подход показал свою эффективность при изучении большого разнообразия явлений, таких как радиационные повреждения, разрушение и развитие биологических структур [60], в частности он позволил смоделировать действие ударной волны на фосфолипиды [61].

Подводный взрыв классифицируют как взрыв на малой глубине (неглубокий), когда глубина воды является малой по сравнению с размерами образованного им кратера. В данном случае процесс образования волны является нелинейным, активно-диссипативным, разительно отличающимся от глубоководного взрыва. На Рис. 1.2 показан процесс распространения волн при взрыве на малой глубине. В его верхней части показан уже заключительный процесс выброса поверхностного слоя воды, с образованием купола брызг и кратера максимального радиуса, а также отходящей от его края диссипативной ведущей волны. Затем данный кратер схлопывается, сдвигая в стороны массы воды и образуется бор (от старо-норвежского bara – волна, зыбь – аномально высокая приливная волна, возникающая в устьях некоторых рек и узких заливов, и движущаяся вверх по их рукавам), за фронтом которой следуют несколько уединенных волн (солитонов).

В эпицентре взрыва, ввиду всех происходящих явлений (подводная ударная волна и потоки воды ломают и разносят объекты на мелкие части), подводный взрыв оказывает разрушительное воздействие на все близлежащие объекты. Образовавшиеся, в результате подводного взрыва, на поверхности воды волны могут долго не затухать и тем самым распространяться на большие расстояния. Это свойство позволяет им воздействовать на объекты, находящиеся на значительРисунок 1.2: Схематическое изображение распространения волн при взрыве на ном удалении, за пределами воздействия подводной ударной волны и вызывать разрушительные цунами. Тем самым волны, порожденные подводным взрывом, тоже могут оказывать разрушительное воздействие на различные объекты, находящиеся на побережье в непосредственной близости от места выхода волны на берег.

1.5 Применение подводных взрывов Данный процесс интересен не только в связи с перспективами и возможностями его непосредственного применения в военно-морской технике, так и с целью его мирного применения. Например, уже сейчас ученые задумываются о применении нано-взрывчатки, которая может помочь в клеточной доставке лекарств и терапии рака [62]. Путем смешивания разных наноматериалов – топлива и окислителя – можно произвести взрыв в нанометровом диапазоне длин, который порождает взрывную волну, распространяющуюся со скоростью до Махов. В состав нано-взрывчатки в качестве топлива входит разреженный порошок наноштырей оксида меди, а в качестве окислителя – наночастицы алюминия. Нано-термит вводится в организм инъекцией и распространяется естественным образом по телу больного, в основном вблизи раковых опухолей. Затем при помощи микрочипа, управляющего взрывом, в место опухоли подается мощный импульс. Ударные волны, сгенерированные этим импульсом, приводят к образованию крошечных проколов в мембранах опухолевых клеток, таким образом открывая их для доставки лекарств. Теми же сверхзвуковыми ударными волнами лекарство может доставляться в клетки опухоли.

Взрывные процессы, как в воде, так и в грунтах нашли применение при ведении дноуглубительных и руслоочистительных работ; строительстве и реконструкции инженерных сооружений (пирсов, причалов, портов, гидростанций);

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

1.6 Математические модели распространения нелинейных волн В современном мире подавляющее большинство явлений являются нелинейными по своей природе; подводный взрыв в данном случае не является исключением. Уравнения Кортевега-де Вриза и уравнения Бусинеска – наиболее распространенные математические модели, которые описывают распространение волн в однородных нелинейных средах со слабой дисперсией. Это моделипортреты волн на мелкой воде. Приближение мелкой воды – это условие малости глубины слоя по сравнению с длинами рассматриваемых волн. Далее в подразделе приведен вывод уравнений Бусинеска и Кортевега-де Вриза, для случая плоских волн в слое воды с плоским дном [63].

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

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

Сложнее вывести граничные условия для поверхности слоя. Уравнение поверхности слоя задается в виде где 0 – глубина слоя в отсутствие волн. Так как частицы воды не могут пересекать поверхность раздела (в силу самого определения этой поверхности), то на поверхности раздела можно положить / = (,, ), / = (,, ).

Тогда, дифференцируя (1.3) по времени, получим Кроме того, если пренебречь поверхностным натяжением воды и считать давление воздуха воды постоянным и равным 0, то на поверхности нужно положить Чтобы записать это условие через скорости и, предположим, что движение воды является потенциальным, т.е.

Тогда уравнения (1.1) можно проинтегрировать и найти :

Из (1.7) и (1.4) следует, что Дифференцируя это равенство по и учитывая (1.6), получим уравнение связывающее компоненты скорости и на поверхности воды:

Учитывая условие потенциальности, уравнения (1.1) можно свести к виду Эти уравнения требуется решить с граничными условиями (1.2), (1.4), (1.9) Если глубина слоя мала по сравнению с длинами рассматриваемых волн (приближение «мелкой» воды), то решение уравнений (1.10) можно искать в виде ряда по координате :

Подставляя (1.11) в (1.10) и приравнивая члены при одинаковых степенях, находим Из граничного условия (1.2) следует, что 0 (, ) 0. Учитывая это, из (1.12) получаем, что (,, ) содержит только четные степени, а (,, ) – нечетные.

Кроме того, как следует из (1.11), Подставляя теперь (1.13) в граничные условия (1.3) и (1.8), получаем В правых частях уравнений (1.15) выписаны члены, имеющие более высокий порядок малости, чем члены в левых частях. Для дальнейшего упрощения уравнений (1.15) предположим, что амплитуды волн достаточно малы, так что в правых частях уравнений (1.15) можно пренебречь нелинейными членами. Тогда получим некий вариант уравнений Бусинеска:

Чаще, однако, уравнениями Бусинеска называют уравнения вида где () – некоторая заданная функция переменной. Уравнения такого вида в одномерном случае можно приближенно получить из (1.15), если в первом уравнении (1.15) положить 0 = 0, а в последнем члене второго уравнения заменить / на /. Тогда получим Чтобы записать уравнения (1.16) в форме уравнений Лагранжа, введем потенциал, так что u =. Тогда уравнения (1.16) можно преобразовать к виду Им соответствует лагранжиан Отсюда находим плотность и поток энергии Полагая в уравнениях (1.16) = 0 +, где 0 – постоянная составляющая переменной, считая величины и u малыми и разлагая функцию () в ряд по : () = (0 ) + 1 + 2 2 +..., можно приближенно переписать уравнения (1.16) в следующем виде:

где 1 и 2 – коэффициенты разложения функции () в ряд по : () = (0 ) + 1 + 2 2 +.... Исключая теперь из этих уравнений переменную u, получим для следующее уравнение:

Путем подбора масштабов зависимой и независимой переменных уравнение (1.22) в одномерном случае может быть перезаписано в виде:

Такое уравнение можно приближено получить и из (1.19), что и было сделано Бусинеском.

Уравнение (1.19) описывает плоские волны, распространяющиеся в обе стороны. Если рассматривать волны, распространяющиеся только в одну сторону, то из уравнений Бусинеска (1.16) следует уравнение Кортевега-де Вриза. Далее продемонстрировано, что это действительно так.

Если пренебречь правыми частями в уравнениях (1.16), то эти уравнения будут иметь решение в виде простых волн. В самом деле, полагая = (0 ), получаем Сравнивая эти уравнения, находим уравнение для (0 ):

решение которого имеет вид Подставляя (1.25) и (1.26) во второе уравнение (1.24), получим для 0 следующее уравнение: ( ) Уравнения (1.16) уже не имеют решения в виде простых волн. Однако, учитывая малость правых частей этих уравнений, из второго уравнения (1.16) в первом приближении можно получить, что:

где (0 ) определяется уравнением (1.26). Подставляя теперь (1.28) в первое уравнение (1.16) и оставляя только члены первого порядка малости, получаем для 0 следующее уравнение:

Перейдя в уравнении (1.29) к новым переменным = (3/2)0, = 0, = и опуская штрихи, получаем уравнение КдВ:

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

Стоит отметить, что уравнение (1.30) впервые было введено Д. Кортевегом и Г. де Вризом в 1895 г. [64]. Как следует из приведенного вывода уравнения КдВ, для волн на «мелкой» воде это уравнение описывает изменение горизонтальной составляющей скорости воды в системе координат, движущейся со скоростью 0. Форма поверхности воды, т.е. форма распространяющейся волны, определяется через уравнения КдВ посредством уравнений (1.26) и (1.28).

1.7 Задача о сильном взрыве Как отмечается в книге Седова [52], при сильном взрыве, которым и является исследуемый в данной работе подводный взрыв, область возмущенного движения отделена от невозмущенного состояния ударной волной. В данном случае при сильном взрыве можно пренебречь давлением перед ударной волной по сравнению с давлением за ударной волной. Учитывая, что 1 = 0, условия на ударной волне записываются в следующем виде:

где – скорость распространения ударной волны, 1 – скорость звука в невозмущенной среде, – показатель адиабаты Чем интенсивнее ударная волна, тем меньше отношение 1 /. На Рис. 1. показана зависимость 1, 2 и 3 от 1 /, согласно которой если в соотношениях (1.31) положить 1 / = 0, 1 = 2 = 3 = 1, то при 1 / < 0.1 можно допустить, что ошибка в значениях 1,2 и 2 менее 5 процентов. Условия на ударной волне при этом примут форму Далее, примем за основные размерные постоянные величины 1 и /1, воспользовавшись указанной выше постановкой задачи о сильном взрыве и уравнеРисунок 1.3: Зависимость величин 1, 2, 3 от соотношения 1 / ниями движения в следующей форме:

где – показатель адиабаты, = 1 для плоских, = 2 для цилиндрических и = 3 для сферических волн. При этом – некоторая постоянная, имеющая ту же размерность, что и энергия 0, выделяющаяся при взрыве; ее размерность равна: [] = 1 2 Видно, что постоянная пропорциональна 0 : 0 =, где – постоянная. Единственный безразмерный параметр в этом случае равен Закон движения ударной волны устанавливается, не решая уравнений движения.

Уравнения движения могут быть различными, однако эти уравнения не должны содержать новых существенных физических постоянных с размерностями, не зависящими от 1 и. В частности, в уравнениях (1.33) нет нужды предполагать, что коэффициент = / имеет постоянное значение.

Для ударной волны координата 2 есть функция времени, а так как из размерных величин, 1 и нельзя образовать безразмерной комбинации, то где * = ; * можно положить равным любому отличному от нуля числу и вычислить значение в зависимости от величины энергии заряда 0. Для простоты можно положить * = 1. Этим на основании решения уравнений движения определится постоянная в формуле 0 =.

Таким образом, для закона движения и для скорости ударной волны получаем в случае сферической симметрии в случае цилиндрической симметрии для полских волн Эти формулы показывают что закон затухания ударной волны зависит от формы заряда. Стоит отметить что формула (1.36) хорошо согласуется с данными взрыва атомной бомбы в Нью-Мексико в 1945 году. Граница ударной волны имела сферическую форму, была резко очерчена и совпадала с ударной волной.

С увеличением времени ударная волна ослабляется, температура за ее фронтом уменьшается. Однако видимость фронта волны при некоторых условиях сохраняется за счет скачка плотности. Опытные данные хорошо подтвердили закон распространения ударной волны (1.36), вывод которого произведен Седовым [52] с помощью теории размерности.

1.8 Молекулярная модель С момента зарождения Максвеллом и Больцманом [65] основ молекулярнойкинетической теории прошло уже более ста тридцати лет. С тех пор достигнут значительный прогресс в обосновании и понимании уравнения Больцмана, создан математический аппарат для описания макроскопической гидрогазодинамики – континуальные уравнения Эйлера и Навье-Стокса, которые применяются в математическом моделировании многих природных явлений и задач [66–68].

Однако, при этом важно отметить, что,например, описание течений вполне удовлетворительно может быть проведено в рамках континуальных уравнений лишь только в областях, где состояние близко к равновесному, т.е. характерные размеры и времена значительно больше кинетических, однако использование данных уравнений не корректно в областях внутри течения с большими градиентами газодинамических параметров. Известно также, что уравнения Навье — Стокса неточно описывают структуру ударной волны из-за сильной неравновесности статистического распределения молекул по энергиям в области резкого изменения параметров газа и д.р. [69, 70]. Это относится также и к описанию кинетики физико-химических процессов во фронте ударной волны. В то же время практические потребности в науке и технике (военной, космической, авиационной, вакуумной, химической и т.д.) требуют решения все более сложных и зачастую нелинейных задач, вне областей применимости классической гидрогазодинамики.

Существуют различные подходы к описанию течений жидкости и газа [71, 72]. Исследуемую задачу можно описывать не только на макроскопическом [73], но и на микроскопическом уровне [74–77]. В отличие от рассмотрения среды как единого континуума, микроскопический молекулярный подход рассматривает корпускулярную структуру среды как миллиарды дискретных молекул и содержит информацию о положениях и скоростях каждой молекулы в любой момент времени. Континуальное описание справедливо до тех пор, пока физически малый объем содержит еще достаточное количество молекул для установления средних величин. И как уже отмечалось выше, уравнения НавьеСтокса не годятся в случае больших градиентов, когда определяемые ими масштабы длины становятся того же порядка, что и среднее расстояние проходимое молекулами между столкновениями(длина свободного пробега). [78, 79]. Отношение длины свободного пробега, к характерному размеру, является числом Кнудсена (1.39).

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

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

Глава Постановка задачи Изучаемая система представляет собой жидкость Леннард-Джонса в двумерном пространстве с периодическими или жесткими границами в продольном направлении и стенкой на нижней границе. На частицы жидкости действует внешнее гравитационное поле в направлении, перпендикулярном поверхности жидкости. Уравнение движения жидкости имеет вид (2.1) где = (0, · ), – гравитационная постоянная. Параметры, использованные в потенциале Леннард-Джонса, соответствуют таковым для воды:

= 0.317, = 0.6502 [80].

2.1 Моделирование подводного взрыва В проблеме моделирования подводного взрыва – его генерации, распространении и взаимодействии волн с мишенью, можно выделить 3 основные части:

Как смоделировать взрывчатку;

Как смоделировать окружение (в нашем случае воду);

Как смоделировать взаимодействие с мишенью;

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

2.1.1 Моделирование окружения В качестве окружения рассматривается вода. Предыдущие компьютерные эксперименты, проведенные с водой [81, 82] показали, что структура ее молекулы сложнее, чем ранее предполагалось. Было найдено, что для правильного описания такого объекта необходимо использовать теорию зарядов связей, где в дополнение к межатомному взаимодействию следует учитывать электрическое взаимодействие между зарядами ковалентных связей [83–85]. Поскольку такой глобальный подход ведет к большим трудностям вычислительного плана и к усложнению модели, было решено построить более простую модель воды, заменив молекулу одной частицей. Взаимодействие таких частиц описывалось с помощью потенциала Леннард-Джонса, рассмотренного подробно в работе [86].

2.1.2 Моделирование ВВ Проблема моделирования взрывчатого вещества, как таковая, представляет собой уже отдельную задачу исследования. Детонация взрывчатых веществ является сложным процессом выделения энергии с образованием ударной волны в реагирующем веществе с вовлечением в этот процесс новых порций взрывчатого вещества, вместе с протеканием химических реакций. Ранее было установлено, что при исследовании плоских ударных волн в воде, взрывчатка может быть исключена [61]. Как результат, получается только две системы - мишень и окрестность вокруг неё. В моделировании подводного взрыва необходимо определить, как будет задаваться взрывчатка. В основу построения модели взрывчатки положен тот же принцип, что и при выборе модели окружения. Вместо компьютерного моделирования на атомарном уровне, сложная модель молекулы взрывчатки заменяется одной частицей. Переход к такой модели значительно упрощает модель. Данный подход используется при моделировании, например, полимерной взрывчатки. При моделировании таких явлений не требуется знать все подробности о строении атомов, достаточно знать, сколько атомов действует как единое целое, формируя кластер.

В данной работе в круговую область заданного размера (с определенными координатами) помещалось известное количество частиц, подвергнутых равномерному сжатию. Вследствие изменения расстояния между частицами данная структура накапливает большое количество упругой энергии (Рис. 2.1). В результате декомпрессии данной области частицы приобретают огромную скорость и, как следствие этого, высокую температуру. Таким образом осуществляется моделирование взрыва. Данная методика нашла отражение и в других наших работах [42, 48, 50, 51].

Рисунок 2.1: Моделирование взрывчатки: ее сложные молекулы заменяются частицами в виде шариков; пунктирными стрелками показано воздействие сил всестороннего сжатия; сплошной стрелкой обозначен переход системы в новое состояние с высвобождением большой энергии 2.2 Расчетная область. Расположение частиц.

Вначале жидкость полностью заполняет нижнюю часть расчетной области ABCD, которая изображена на Рис. 2.2.

Размеры ABCD: 0 < 1 <, 0 < 2 < в продольном и поперечном направлениях. Частицы расположены упорядоченно на расстоянии 2 6 друг от друга, что соответствует минимуму энергии Рис. 2.3.

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

Рисунок 2.4: Расположение частиц в кольцевом заряде Здесь 1, 2 – координаты центра заряда, – расстояние в радиальном направлении, на котором располагаются кольца, – полярный угол между частицами одного из внутренних колец. Далее проводится декомпрессия частиц в заряде, имитирующая подводный взрыв.

2.3 Граничные условия В данной постановке, на вертикальных границах AB, CD (см. Рис. 2.2) задавались как периодические граничные условия, так и граничные условия, соответствующие твердой границе. При компьютерном моделировании периодические граничные условия часто применяются при моделировании большой системы, путем использования небольшой её части.

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

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

При использовании потенциала Леннарда-Джонса, дополнительная репульсивная сила действующая на частицу выражается формулой (2.2) Здесь обозначает расстояние между частицей и стенкой. Данная сила имеет только репульсивную составляющую, так как действует на частицу при расстояРисунок 2.6: Граничное условие, реализующее твердую стенку нии меньшем 2 6, данное расстояние соответствует минимуму энергии. Следует отметить, что вектор расстояния между частицей и её отображением : = = (0, 2), что вносит отрицательный знак. Необходимо также учитывать, что временной шаг должен быть выбран достаточно малым, для того чтобы гарантировать стабильность метода.

2.4 Методика наблюдения динамики ударной волны Немаловажной особенностью метода молекулярной динамики является изучение одновременно как динамической структуры системы, так и наблюдение за движением отдельных её частиц, вместе с расчетом усредненных характеристик системы: энергии, давления, температуры и др., особенно с учетом их распределения и изменения в системе с течением времени. Специально для наблюдения за динамикой зарождения и распространения ударной волны была разработана и применена специальная методика. Идея заключается в том, чтобы через определенные моменты времени рассчитывать радиально распределенную плотность частиц вдоль линии направления движения ударной волны, в заранее определенной области пространства, по нескольким значимым направлениям. Таким образом в выходном файле оказываются данные зависимости распределения радиальной плотности от координат в различные моменты времени. После построения графиков на основе этих зависимостей возможен дальнейший анализ динамики данного процесса, определения скорости движения ударной волны, ее интенсивности и затухания, зависящих от формы и мощности моделируемого заряда. На Рис. 2.7 схематично представлено разбиение расчетной области на 3 сектора по 30 градусов в диапазоне от 0 до 90. В расчетную область, в центр заряда помещается начало системы координат. Сектор данной системы координат, соответствующий положительным направлениям координатных осей, также дополнительно разбивается на три подсектора по 30 градусов каждый. В Рисунок 2.7: Схематическое изображение разбиения области распространения радиальном направлении задается шаг, отвечающий за вторую составляющую разбиения области. В итоге мы имеем систему подсекторов окружностей, в каждом из которых, через заранее определенный шаг по времени будет производиться расчет радиальной плотности частиц по следующей формуле:

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

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

2.5 Алгоритм решения задачи Упрощенно, алгоритм решения задачи можно описать следующим образом:

Производится ввод данных, таких как: температура системы, параметры потенциала и, радиус его обрезания, величина силы тяжести, а также масса, число частиц и их координаты, шаг по времени, время окончания расчета ;

Начальные скорости частиц задаются, следуя распределению МаксвеллаБольцмана;

Запускается вычислительный процесс, в котором по известным координатам частиц системы в начальный момент времени, вычисляются силы, которые действуют на каждую частицу По известным силам вычисляются новые координаты, силы и скорости на новом временном шаге + ;

Далее процесс повторяется до момента достижения времени окончания расчета ;

Блок-схема вычислительной программы, при условии выбора метода НордсикаГира для решения дифференциального уравнения, представлена ниже:

Рисунок 2.8: Блок-схема вычислительной программы 2.6 Цели и задачи.

В качестве метода исследования было выбрано численное моделирование.

При этом были поставлены следующие задачи:

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

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

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

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

Рассмотреть динамику подводного взрыва;

Провести анализ процесса распространения ударных волн;

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

Глава Вычислительные аспекты 3.1 Метод молекулярной динамики Метод молекулярной динамики – один из методов компьютерного моделирования, использующий классическую механику для изучения поведения системы с течением времени. При использовании данного метода численно решается система уравнений (3.1) Интегрирование второго закона Ньютона с помощью различных схем, от простейших методов Эйлера первого порядка до схем предиктор-корректор, позволяет определить траектории атомов. Наиболее важными требованиями для таких схем является устойчивость, точность и высокая скорость работы. Сила, действующая на атом, определяется как отрицательная производная потенциальной энергии (3.2) Если заданы начальные координаты и скорости частиц, то эволюция системы во времени зависит только от потенциала, который определяет взаимодействие между атомами.

3.2 Потенциал Леннард-Джонса Потенциал Леннард–Джонса [87–89] (потенциал 6-12) выражает простую модель парного взаимодействия неполярных молекул и описывает зависимость энергии взаимодействия двух частиц от расстояния между ними. Он выражается формулой (3.3). Такой потенциал хорошо описывает взаимодействие сферических неполярных молекул и поэтому широко используется в расчетах, а также при компьютерном моделировании здесь – расстояние между центрами частиц и, – глубина потенциальной ямы, – расстояние, на котором энергия взаимодействия становится равной нулю. Параметры и являются характеристиками вещества. На Рис. 3.1 показан общий вид потенциала Леннарда-Джонса при параметрах = 1 и = 1.

При больших значениях молекулы притягиваются, что соответствует члену в формуле (3.3). Эту зависимость можно обосновать теоретически и показать, что обусловлена она силами Ван-дер-Ваальса (диполь-дипольное индуцированное взаимодействие). На малых же расстояниях молекулы отталкиваются из-за обменного взаимодействия (при перекрытии электронных облаков) Потенциальная функция для частиц приобретает вид двойной суммы Соответственно сила, действующая на частицу, определяется (3.5) как градиент по направлению здесь = направление вектора между частицами и в точках и. Результирующая сила (3.6), действующая на частицу, состоит из суммы сил между частицей и частицами Рассматриваемый потенциал действует на бесконечном расстоянии между частицами, что представляет сложность при численном моделировании, так как необходимо учитывать взаимодействие между всеми частицами, что приводит к резкому росту времени счета при увеличении числа частиц. Такой потенциал также создает проблемы при использовании периодических граничных условий, так как при малой длине стороны моделируемого объема частица взаимодействует сама с собой. Как видно из Рис. 3.1, потенциал (а с ним и действующая сила) спадает очень быстро с увеличением расстояния между частицами. Именно по этим причинам на практике используется усеченный потенциал.

Мы пренебрегаем взаимодействием частицы с частицей, если расстояние оказывается больше определенного порогового значения. Потенциал ЛеннардДжонса (3.3) при этом аппроксимируется выражением:

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

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

Тем не менее, зачастую на практике для описания взаимодействия между частицами используется модифицированный потенциал Леннард-Джонса вида:

где ( ) сглаживающая функция, определяемая следующим образом:

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

Предположим теперь, что частицы более или менее равномерно распределены по всей расчетной области. Тогда может быть выбран, чтобы количество элементов, оставшихся в усеченной сумме в (3.8) и (3.9), ограничено не зависело от числа частиц. Сложность вычисления потенциала и сил теперь пропорциональна, то есть порядка ( ). Это приводит к значительному сокращению вычислительных затрат по сравнению со сложностью порядка ( 2 ), которая будет в случае стандартного подхода вычисления силы по всем парам частиц.

Рисунок 3.2: Потенциал Леннард-Джонса со сглаживающей функцией и без 3.3 Дискретизация уравнений движения.

3.3.1 Метод Верле Для получения траекторий атомов в пространстве и времени существуют различные численные схемы решения системы классических уравнений движения, от простейших методов Эйлера первого порядка, до схем предикторкорректор высоких порядков точности [90–92]. Одним из наиболее распространенных и в тоже время стабильным и эффективным подходом к временной дискретизации уравнений Ньютона является алгоритм Верле [93,94]. Данный метод заключается в следующем. Мы разбиваем временной интервал [0, ], на котором система дифференциальных уравнений (3.1) должна быть решена, на подынтервалов одинакового размера, = /. В этом случае мы получаем сетку, состоящую из точек = ·, = 0..., находящихся на концах подынтервалов. Дифференциальное уравнение теперь рассматривается в этих временных точках. В системе обыкновенных дифференциальных уравнений второго порядка (3.1), дифференциальное отношение может быть заменено разностным отношением в каждый момент времени, = 0,..., 1. Вторая производная 2 /2 аппроксимируется в точке следующим дифференциальным оператором: [ 2 ] С учетом переобозначений = ( ), = ( ), = ( ) получаем:

Скорости, как производные по координате, могут быть аппроксимированы центрально-разностной схемой:

В нашей программе использовалась видоизмененный вариант метода Верле второго порядка точности, т.е. (2 ). Если мы выразим из (3.14) 1 и подставим результат в (3.13), тогда получим выражения для Кроме того, (3.13) и (3.15) дают Добавляя соответствующее выражение для, получаем:

И далее, используя (3.14), Уравнение (3.15) вместе с уравнением (3.18) называется скоростным методом Верле (Velocity Verlet method).

На Рис. 3.3 схематично представлена последовательность шагов при временной итерации. Таким образом, зная начальное распределение координат и скоростей, = 1,...,, в начальный момент времени = 0, мы на каждом временном шаге до момента окончания вычислений последовательно вычисляем необходимые нам значения согласно (3.15), (3.9) и (3.18).

3.3.2 Метод Нордсика-Гира При написании современных научных вычислительных программ всегда остро стоит вопрос об интегрировании системы обыкновенных дифференциальных уравнений, выбор метода в данном случае играет одну из ключевых ролей. Располагая начальной конфигурацией системы (позиции всех атомов составляющих данную систему) в начальный момент времени = 0, мы пытаемся вычислить ее конфигурацию в момент времени +. В итоге, вычисляя ее состояние на каждой новой итерации, мы имеем представление об эволюции данной системы. Точность решения при этом строго зависит от примененного нами вычислительного алгоритма и размера временного шага. Метод решения системы ОДУ должен в идеале сочетать в себе 2 важнейших составляющих – это скорость работы и обеспечение должного порядка точности. Как правило, выбор метода всегда является своего рода компромиссом. В данной задаче наряду с методом Верле, был использован также метод Гира [95,96] в представлении вектором Нордсика [97] 6 - ого порядка точности. Преимущество данного метода заключается в его надежность и высоком порядке точности(до 6-ого порядка) при относительно умеренных(больших) временных шагах.

Ветор Нордсика представляется следующим соотношением:

В данном случае означает аппроксимацию - ого порядка в точке, где () является точным решением дифференциального уравнения Когда рассматривается дифференциальное уравнение вида (3.20) и вектор задан, значения +1 вычисляются, исходя из разложения его компонент в ряд Тейлора(пример для = 3):

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

предикторный шаг в матричной форме записи будет иметь вид, где r являются следующим образом смасштабированными производными по времени: r1 = (r0 /), r2 = (1/2)2 (2 r0 /2 ), r3 = (1/6)3 (3 r0 /3 ), r4 = (1/24)4 (4 r0 /4 ). Матричная форма записи делает предикторный шаг легко применимым. При этом число коэффициентов в разложении Тейлора определяет порядок метода. Разница между значениями второй смаштабированной производной по времени r на предикторном шаге и значением r, полученным путем подставления значения r0 в уравнения движения, определяется формой:

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

Коэффициенты являются переменными. Их значения влияют, как на устойчивость численного алгоритма, так и на точность полученных результатов. Гир в работах 1966, 1971 годов описал и предложил оптимальные коэффициенты для дифференциального уравнения (3.22) второго порядка: 0 = 3/20, 1 = 251/360, 2 = 1, 3 = 11/18, 4 = 1/6, 5 = 1/60.

Как и в случае с методом Верле, временной интервал [0, ], на котором система дифференциальных уравнений (3.1) должна быть решена, разбивается на подынтервалов одинакового размера, = /. Получается одномерная сетка, состоящая из точек = ·, = 0..., находящихся на концах подынтервалов. Дифференциальное уравнение теперь рассматривается в этих временных точках. С учетом обозначений = ( ), = ( ), = ( ) и вектора Нордсика вида:

получаем следующий систему уравнений на предикторном этапе:

На корректорном этапе:

Важно отметить, что в реализованной программе значения для скоростей вычисляются уже в явной форме согласно представленной выше системе уравнений (3.27), (3.28), а не в виде второго члена в векторе Нордсика, которым является второй член в разложении в ряд Тейлора. Данной системой уравнений, в два этапа, на каждой итерации определяются новые значения переменных в “модифицированном“ векторе Нордсика. Значения коэффициентов 0,..., приведены согласно Гиру [95].

3.4 Оптимизация вычислений 3.4.1 Метод связных ячеек Организация вычислительных данных для данной частицы, так чтобы её взаимодействие с соседними частицами было найдено наиболее эффективным методом, оказывается актуальным при исследовании систем, состоящих из нескольких десятков тысяч частиц. При моделировании систем такого размера мы ставим задачу добиться оптимизации вычислений, с целью уменьшения компьютерных затрат. В данной работе для этого использовался метод связанных ячеек (Linked Cell Method), который сравнительно легко осуществим и в то же время очень эффективен [98–100]. Данный метод применяется для вычисления сил и энергий быстро убывающих потенциалов, в случае которым как раз и является потенциал Леннард-Джонса.

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

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

Суммы в (3.13) или (3.9) теперь разбиваются на частичные суммы, соответствующие декомпозиции расчетной области на ячейки. Теперь сила действующая на частицу в ячейке приобретает вид:

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

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

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

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

Рисунок 3.5: Метод связных ячеек. Расчетная область с периодическими Сложность вычисления сил, действующих на частицу составляет · операций, если частицы распределены почти равномерно. По сравнению с вычислением силы по всем парам частиц, сложность уменьшается c ( 2 ) до ( ).

Как утверждалось выше, ячейки должны быть, как минимум, радиуса обрезания потенциала, для возможности применения данного метода. В данной работе мы определяли продольный и поперечный размер прямоугольной области как - 1, 2 и также число ячеек вдоль и поперек - 1, 2. Из соотношения /, = 1, 2 следует, что наибольшее число ячеек в продольном и поперечном направлении может быть вычислено как:

Здесь, [] определяет наибольшее целое. 1, 2 – продольное и поперечное направление.

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

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

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

Далее используем длину, энергию и массу для масштабирования переменных следующим образом:

и далее путем подстановки в (3.33):

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

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

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

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

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

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

Значение энергии системы, состоящей из частиц и находящейся при температуре, определяется формулой (3.48):

где = 1.38065 · 1023 Дж · 1 – постоянная Больцмана. Мы можем привести Таким образом:

В соответствие с этим, имеет место формула (3.51):

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

3.7 Точность метода молекулярной динамики При применении того или иного численного метода для решения задачи компьютерного моделирования принципиально стоит вопрос о его точности.

Он приобретает первостепенное значение, поскольку результаты моделирования приобретают силу экспериментального факта в тех случаях, когда моделируются системы, экспериментальное изучение которых затруднено или попросту невозможно. Данный вопрос подробно изучался в работах Нормана [103, 104].

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

Ошибки округления;

Ошибки, связанные с использованием различных схем и алгоритмов численного интегрирования уравнений Ньютона;

Ошибки, обусловленные конечным размером расчетной области и использованием периодических граничных условий;

Корреляции, обусловленные конечным числом частиц в моделирующей системе;

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

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

Поэтому основным критерием использования той или иной конечно-разностной схемы является прежде всего правильное определение доверительного временного интервала, в течение которого эти ошибки остаются достаточно малыми, чтобы существенно не влиять на качество расчета. Если это время равно, то ошибки определения динамических переменных на таких временах должны быть малыми. Поскольку при использовании конечно-разностных схем того или иного порядка ошибка вычисления динамических переменных пропорциональна шагу дискретизации, то этот шаг должен быть таким, чтобы +1 0. Из условия экстремума (; ) = 0 получаем У данного уравнения существуют два корня таким образом в зависимости от скорости потенциальная энергия имеет три формы (Рис. 5.2). Она равна нулю при = 0 и = 2 1 2. Полный анализ Рисунок 5.2: Влияние скорости распространения ударной волны на эффективную потенциальную энергию (верхняя часть) и на форму безразмерной линейной плотности деформации межатомных связей (нижняя движения для этих трех случаев приведен в [106]. Мы заинтересованы только в случае, когда >, = 0. В данном случае можно записать решение в виде [106].

Здесь = 3( /)2 1, = /, 0 = /. Это решение описывает уединенную волну сжатия или солитон сжатия (Рис. 5.3). Для удобства интерпретаРисунок 5.3: Солитон сжатия ции результатов с точки зрения континуума, перейдем от деформации к перемещению. Поскольку = /, после интегрирования по, мы получаем:

Это решение с точки зрения перемещения имеет форму ударной волны которая движется в положительном направлении оси с постоянной скоростью V. В точке = + inf нет никакого возмущения, поскольку у волны нет времени чтобы дойти до нее, поэтому = 0 и = | | = 0 /. И следовательно где = / ширина ударной волны. Важно отметить, что 0 2 = =, т.е. произведение амплитуды ударной волны на квадрат ее ширины является константой. С уменьшением скорости волны, ее амплитуда также уменьшается, поэтому в пределе волна оказывается полностью смазаной. Стоит отметить, что скорость связана с амплитудой 0 соотношением (5.38) Так как = 2 (0), мы можем сказать, что физический смысл пространственной дисперсии волн связан с кубической ангармоничностью потенциала межатомного взаимодействия. Для справки следует отметить, что уравнение КдВ на атомарном уровне было рассмотрено для конкретных форм межатомного потенциала, а именно для цепочки Тода [114].

5.2.1 Уравнение Кортевега-де Вриза и уравнение Шредингера Здесь мы будем следовать подходу изложенному в [111]. Пусть функция (, ) удовлетворяет уравнению Шредингера:

где H - оператор Гамильтона. Квантовая механика может быть сформулирована и другим способом, при котором зависимость от времени с волновых функций перенесена на операторы (представление Гейзенберга). Если – Эрмитов оператор, удовлетворяющий некоторой физической величине, тогда производная этого оператора по времени определяется соотношением (5.40) где [, ] = коммутатор. Если оператор не зависит от времени, то все его собственные числа, которые определяются соотношением (5.42) также не зависят от времени. Уравнение / = 0 является условием постоянства оператора и, следовательно, Операторы и определяются следующим образом записан в форме (5.44) В результате прямых и довольно кропотливых вычислений получаем:

Перепишем уравнение = в форме Мы можем интерпретировать это уравнение как одномерное стационарное уравнение Шредингера, где первый член является кинетической энергией, является полной энергией, а потенциальной энергией. В математике данное уравнение известно как уравнение на собственные значения [114]. Для функции (, 0), которая быстро уменьшается при ||, мы имеем конечное число дискретных, отрицательных собственных значений [115]. В качестве примера, рассмотрим начальное значение из [114] и рассмотрим временную эволюцию решения уравнения КдВ. Когда = 1, и амплитуда равна 2, форма волны не меняется со временем и она движется вправо.

Это солитон или устойчивая уединенная волна, которую наблюдал и исследовал Скотт Расселл. При = 2 и амплитуде 6, собственные числа = 0, 1, 4. В результате, исходное состояние распадается на два солитона, один с собственными значениями = 0, 1, и другой, имеющий = 0, 4. Вся энергия и импульс начального состояния уносятся этими двумя солитонами. Когда не является целым числом, у нас получаются солитоны и некоторая рябь (осциллирующий хвост). Они возникают из начального состояния и движутся вправо, в то время как вся рябь распространяется влево с убывающей амплитудой.

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

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

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

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

Здесь – скорость уединенной волны, – ускорение свободного падения, – глубина воды, 0 – максимальный подъем жидкости, над уровнем воды. Среди многих исследований, Кортевег и де Вриз рассматривали движение мелких волн на воде. Они нашли отражение в уравнении (5.51) где – подъем уровня водной поверхности, =, – поверхностное натяжение, – плотность воды. Записывая уравнение КдВ в канонической форме мы сводим двумерное явление к одномерному. Так как () пропорционально дополнительной массе, можно рассматривать как изменение линейной плотности среды. Рассматривая уравнение КдВ на атомарном уровне, мы получили следующее уравнение для стоячей волны где – безразмерная линейная плотность деформации межатомной связи. Поскольку деформация межатомной связи ведет либо к увеличению, либо к уменьшению межатомного расстояния, можно рассматривать () как безразмерное изменение линейной плотности среды. Было показано, что постоянная Грюнайзена = 2 (0). Поскольку 1/, можно придать физический смысл дисперсии, связав ее с кубической ангармоничностью потенциала межатомного взаимодействия. Осталось обсудить возможную структуру функции (, ). Согласно первоначальной структуре уравнения КдВ, она имеет два параметра, которые могут изменяться в зависимости от условий эксперимента, плотности и глубины. В свою очередь эти обе величины связаны с давлением. Если ударная волна распространяется, когда 0, плотность и глубина почти не меняются и поэтому коэффициент почти постоянен, поэтому и форма ударной волны почти постоянна, по крайней мере, после некоторой начальной стадии. Такая ситуация наблюдается для ударных волн распространяющихся в секторе 0 30.

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

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

5.3 Выводы Был разработан метод изучения структуры ударных волн, возникающих при подводных взрывах и наблюдаемых в рамках молекулярно-динамических расчетов. Наиболее представительным результатом исследования является то, что форма ударной волны меняется со временем, как если бы она распадалась на две части. Данное явление было проанализировано при помощи подхода, применяемого для явления солитонов. Получено уравнение для двумерных ударных волн, обобщающее уравнение КдВ для одномерной среды. В полярных координатах оно имеет вид + + (, ) = 0. Показано, что дисперсия c кубическим ангармонизмом потенциала межатомного взаимодействия.

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