Редактирование: Совершенствование алгоритмов численного моделирования в методе динамики частиц

== Задача ==

* ...

'''Цели''':
* Создание программы на базе существующего кафедрального программного пакета по молекулярно-динамическому расчёту для
** Молекулярно-динамическое моделирования кипящего слоя
** Моделирования провала на земной поверхности при перепуске налегающих пород в карстовую полость
** Многого другого...
* Анализ и попытка улучшения эффективности работы и возможностей программы посредством добавления новых численных методов и приёмов счёта.

'''Замечание''':
* Ниже приведены приёмы и методы, уже запрограммированные и активно отлаживающиеся в настоящее время.

== Построение модифицированного метода Рунге-Кутты 4 порядка ==

Рассмотрим задачу Коши

<math> dx/dt=f(x,t) \ \ (1) </math>

<math> x(0)=x_0 </math>

Для неизвестной вектор-функции x(t), в качестве которой для примера может быть взят вектор <math> (r,v) = (x,y,z,u,v,w) </math> координат позиции и скорости тела. Данная задача может быть решена численно классическим методом Рунге-Кутты четвёртого порядка. 

<math> k_1=\Delta t f(x_n,t_n ) </math>

<math> k_2=\Delta t f(x_n+k_1/2,t_n+Δt/2) </math>

<math> k_3=\Delta t f(x_n+k_2/2,t_n+Δt/2) </math>

<math> k_4=\Delta t f(x_n+k_3,t_n+Δt) </math>

<math> x_{n+1}=x_n+1/6 (k_1+2k_2+2k_3+k_4 )  \ \ (2) </math>

По сравнению с методами Эйлера, Лагранжа и Верле, данный метод имеет более высокий порядок точности. Однако классический метод Рунге-Кутты четвёртого порядка имеет одну особенность, связанную с необходимостью вычислять функцию <math> f(x,t) </math> четыре раза за одну временную итерацию. Потому этот метод становится неэффективным в вычислительных задачах, где основное расчётное время тратится на вычисление правой части системы дифференциальных уравнений, как, например, это имеет место в случае расчёта молекулярно-динамической задачи множества частиц. Вследствие данной особенности применение метода Рунге-Кутты становится неэффективным и даже его исключительная точность теряет свою значимость.

Ниже приводится модификация метода Рунге-Кутты 4 порядка, где с помощью одного хитрого приёма удаётся избежать многократного вычисления функцию <math> f(x,t) </math> на одном временном шаге и в то же время сохранить высокий порядок по времени.

Идея заключается в разложении функций <math> f(x_n+k_i/2,t_n+Δt/2) </math> в ряд Тейлора в окрестности точки <math> (x_n,t_n) </math>. 

<math> f(x_n+k_i/2,t_n+Δt/2)= \frac {\partial f} {\partial x} (x_n,t_n )∙k_i/2+  \frac {\partial f} {\partial t} (x_n,t_n )  Δt/2 + ...  \ \ (3)	</math>

Здесь присутствуют малоприятные производные, однако, как потом окажется, с ними можно будет легко разобраться. Сколько членов в разложении нужно оставить, чтобы в схеме сохранился четвёртый порядок? – До <math> (Δt)^4 </math> и <math> (k_i )^4 </math> или меньше?

Для слагаемых с локальными производными по времени ответ очевиден – необходимо удерживать всё вплоть до <math> (Δt)^4 </math>, ибо в противном случае мы потеряем наш 4-й порядок по времени для схемы в целом. Однако для <math> k_i </math> на самом деле достаточно только первой производной. 

В случае, когда правая часть (1) не зависит явно от времени, (3) предельно упрощается.

<math> f(x_n+k_i/2,t_n+Δt/2)=∂f/∂x (x_n,t_n )∙k_i/2 \ \ (4) </math>

Данная ситуация имеет место при молекулярно-динамическом моделировании, поскольку потенциал взаимодействия, как правило, является функцией только координат и скоростей частиц.

Запишем (4) в случае молекулярно-динамического моделирования. В нашем случае неизвестная вектор-функция  <math> x = {r,v} </math>.

<math> (x,t)=(v,\frac{F(r)}{m}) </math>

<math> f(x_n+k_i/2,t_n+Δt/2) = (v,\frac{F(r+k_i^r/2)}{m}) = ( v+\frac {k_1^v}{2},\frac {F(r+\frac {k_1^r}{2})}{m} )= </math>

<math> =( v+\frac {k_1^v}{2},\frac {F_n + \frac {dF}{dr} (r)  \frac{k_1^r}{2}}{m} )  \ \ (5) </math>

где

<math> \frac {dF}{dr} = \frac {\Delta F}{\Delta r} = \frac {\Delta F}{\Delta r} \frac {\Delta t}{\Delta t} = \frac {\Delta F}{\Delta t} \frac {\Delta t}{\Delta r} =\frac {F_n - F_{n-1}}{\Delta t } \frac {1}{v} \ \ (6) </math>

где <math>F_n</math> и <math>F_{n-1}</math> - силы, действующие на рассматриваемую частицу на данном временном слое  <math>n</math> и на предыдущем временном слое  <math>(n-1)</math> соответственно. Это значит, что в памяти необходимо хранить силы с прошлого временного слоя, что не очень хорошо. Это главный, но не очень существенный минус данного метода.
С учётом вышесказанного, мы можем записать конечно-разностную схему Рунге-Кутта для молекулярно-динамического моделирования

<math>k_1^v=\Delta t \frac {F_n}{m}</math>

<math>k_2^v=\Delta t \frac {F_n}{m}+\frac{F_n-F_(n-1)}{2m}</math>

<math>k_3^v=\Delta t \frac {F_n}{m}+\frac{F_n-F_(n-1)}{2m}</math>

<math>k_4^v=\Delta t \frac {F_n}{m}+\frac{F_n-F_(n-1)}{m}</math>

<math>k_1^r=\Delta t \frac {F_n}{m}</math>

<math>k_2^r=\Delta tv_n+\frac {v_n-v_{n-1}}{2}=\Delta t ( v_n+\frac{k_1^v}{2}) </math>

<math>k_3^r=\Delta tv_n+\frac {v_n-v_{n-1}}{2}= \Delta t ( v_n+\frac{k_2^v}{2})</math>

<math>k_4^r=\Delta tv_n+v_n-v_{n-1}= \Delta t(v_n+k_3^v )</math>

<math>r_{n+1}=r_n+1/6 (k_1^r+2k_2^r+2k_3^r+k_4^r )</math>

<math>v_{n+1}=v_n+1/6 (k_1^v+2k_2^v+2k_3^v+k_4^v )	\ \ (8) </math>

Выражения для <math> k_i^r</math> представлены в двух видах: один – с учётом рассмотренного упрощения, а второй - классический Рунге-Куттовский, который применим для уравнения по <math> r </math> ввиду простой правой части.

== Обезразмеривание системы как способ уменьшения накопления вычислительной ошибки ==

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

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

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

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

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

Рассмотрим задачу Коши

<math> \frac{dv_i}{dt}=\frac{f_i (r_{ij},v_{ij} )}{m_i}  </math>

<math> \frac{dr_i}{dt}=v_i </math>

<math> r_i (0)=r_i^0 </math>

<math>v_i (0)=v_i^0 </math>

<math>i=1,..,N</math>

Пусть <math>\tau</math> – характерное время протекания рассматриваемых процессов, <math>a</math> – характерная длина расчётной области, <math>m</math> – характерная масса частиц. 

Сделаем замену переменных

<math>r_i = a R_i</math>

<math> v_i = \frac{a}{\tau} V_i</math>

<math>t = \tau </math>

<math> f_i = m \frac{a}{\tau^2} F_i </math>

<math>i=1,..,N</math>

где <math>\tau, R_i,V_i,F_i </math> – новые безразмерные параметры (время, радиус-вектор, вектор скорости и сила соответственно).

Исходная система Коши в безразмерном виде будет таковой

<math>\frac{dV_i}{d\tau}=F_i (R_{ij},V_{ij} )</math>

<math>\frac{dR_i}{d\tau}=V_i</math>

<math> R_i (0)=(r_i^0)/a</math>

<math>V_i (0)=(v_i^0 \tau)/a</math>

<math>i=1,..,N</math>

Однако на этом обезразмеривание ещё не заканчивается. 

<math> f_i (r_{ij} )=\sum_{j=0}^{N}f(r_{ij})</math>

где

<math> f(r_{ij} )=- \frac{\partial u}{\partial r_{ij}} </math>

Если в качестве потенциала взаимодействия U взять потенциал Леннарда-Джонса

<math> u(r_{ij} )=d[(\frac{a}{r_{ij}} )^{12}-2(\frac{a}{r_{ij}} )^6 ]=d[(\frac{1}{R_{ij}} )^{12}-2( \frac{1}{R_{ij}} )^6 ] </math>

Глубина <math>d</math> потенциальной ямы имеет размерность энергии. Если мы хотим обезразмерить потенциал, то нам достаточно ввести новую безразмерную глубину <math> D </math> потенциальной ямы по правилу

<math>d=D \frac{m a^2}{τ^2}</math>

Тогда новая обезразмеренная потенциальная энергия запишется в виде

<math>U(R_ij )=D ((\frac{1}{R_{ij}} )^{12}-2(\frac{1}{R_{ij}} )^6 )	</math>

А безразмерная сила <math>F</math> будет выражаться  через потенциал привычной формулой

<math>F(R_{ij} ) = - \frac{\partial U}{\partial R_{ij}}	</math>

== Frozen Particles & Press Particles ==

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

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

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

По аналогии с Frozen Layer может быть построен и слой частиц для моделирования "граничных условий" в напряжениях.