10.3. АЛГОРИТМ С СОХРАНЕНИЕМ ЭНЕРГИИ

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

где - объем ячейки, причем это определение справедливо при любом законе взаимодействия, а не только при законе Кулона [Reitz and Milford, 1960]. В частности, равенство (10.7) применимо, когда кулоновская сила сглажена на коротких длинах волн. Плотность заряда в точках сетки определяется как обычно:

Если мы рассчитываем силу, действующую на частицу, по формуле

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

Из (10.7) и (10.9) имеем

Из определения (10.8) получаем для входящей в первое слагаемое, и для расчета требуется только несколько значении Однако при всех Так как это слагаемое дает вклад во все ячейки, его вычисление было бы очень дорогостоящим. Полезно равенство (10.10) переписать в виде

Ниже в этом параграфе показано, что вторая сумма всегда равна нулю. Тогда сила, действующая на частицу, получается из первой суммы выражения (10.11):

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

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

Силу, действующую на частицу, можно затем получить из градиента потенциального поля, как в последнем равенстве (10.12), в котором считается постоянным при дифференцировании Способ, рассмотренный в работе [Lewis, 1970 а], аналогичен нашему, см. § 10.5.

Заметим, что одна и та же интерполяционная функция используется в (10.8) и в (10.12). Левис привел примеры использования линейной интерполяции первого порядка в случае одного и двух измерений. Однако в работах о, 1973] не было никаких ограничений на эти веса. Интерполяция нулевого порядка NGP, в которой является разрывной функцией, не подходит из-за несуществования градиента в соотношении (10.12).

Выше отмечалось, что вторая сумма в (10.11) исчезает; это действительно так, если имеет место (см. задачу 10.2) равенство

где индексы (1) и (2) относятся к двум различным распределениям плотности и соответствующим им потенциалам: В пределе это утверждение переходит в теорему взаимности Грина для электростатики [Jackson, 1975]. Взаимность имеет место, когда является решением разностного уравнения вида

с

(см. задачу 10.3). Симметрия обычно возникает естественным образом при записи уравнения Пуассона в общей криволинейной системе координат. В нейтральной плазме с периодическими граничными условиями и для уравнения Пуассона, симметричного относительного отражения в плоскостях решетки, эта симметрия обеспечена. Формально Первое равенство следует из отражения, а второе — из трансляции на величину В способе Левиса А симметрично для разностного уравнения Пуассона (10.30), так как интеграл инвариантен при перестановке индексов Это последнее свойство вносит существенно меньше ограничений, чем (10.30), и, следовательно, закон сохранения энергии справедлив для существенно более широкого класса алгоритмов, чем это было указано Левисом.

В криволинейных координатах выполняются также условия при а

так как если однородно в пространстве. Эти условия влияют на знак энергии поля (см. задачу 10.5).

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

Можно сделать выводы: когда имеет место взаимность, как в (10.14), сила, используемая в коде с сохранением энергии, идентична градиенту со знаком минус от полной энергии поля. Анализ соотношения (10.23) показывает, что взаимность является необходимым условием.

Используя однородное сплайновое взвешивание порядка (см. § 8.8), из равенства (10.12) получаем выражение для силы

Учитывая тождество для производной

которое следует из определения как огибающей при соответствующем методу NGP, запишем силу в виде

где и (см. задачу 10.6)

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

Задачи

(см. скан)

(см. скан)