LINUX.ORG.RU

вопрос про интегрирование уравнений движения

 , ,


0

2

$Cast AIv

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

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

★★★★★

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

dmfd
()

Я правильно понял: у тебя ODES из двух дифференциальных уравнений. Делаешь всё leapfrog-ом, всё ОК. Делаешь одну часть уравнения аналитически, решение не стабильное?

anonymous
()
Ответ на: комментарий от anonymous

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

1)для каждой частицы решаю уравнения движения за время dt0, интерполируя электрическое и магнитное поле по ближайшему узлу (NearestGridPoint). За время dt0 ЭМ поля считаются постоянными.

2)по необходимости, перебрасываю частицу в соседнюю ячейку

3)интерполирую плотность заряда по узлам, учитывая время нахождения частицы в узле, и решаю ур.е. пуассона и для магн поля.

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

ну так вот если решаю уравнения движения в виде

vx1=-(1.*(-4.*k*Wc*dt0*vy0-4.*vx0+k**2*Wc**2*dt0**2*vx0-4.*k*mu*Ex_*dt0))/(4.+k**2*Wc**2*dt0**2)

vy1=-(1.*(vy0*k**2*Wc**2*dt0**2-4.*vy0+4.*k*Wc*dt0*vx0+2.*k**2*Wc*dt0**2*mu*Ex_))/(4.+k**2*Wc**2*dt0**2)

x1=vx1*dt0+x0

где Wc -- циклотронная частота, k -- знак заряда, в данном случае -1, mu -- отношение заряда к массе, Ex_, Hz_ -- электрическое и магн. поле в точке x0, dt0 -- шаг по времени, который естественно меньше чем циклотронный или плазменный период.

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

однако, можно записать уравнения движения в виде

C0=(mu*Ex(ix1)/Wc+vy0)
C1 = C0*k/Wc
C2 = vx0/Wc
					vx1=k*C0*sin(Wc*dt1)+vx0*cos(Wc*dt1)
xy1=C0*cos(Wc*dt)-k*vx0*sin(Wc*dt)-Ex(ix1)/Hzc(ix1)
x1=x0+C1*(1-cos(Wc*dt1))+C2*sin(Wc*dt1)

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

dt1=(abs(2*atan(C2/C1))+abs(acos((C1)/sqrt(C1**2+C2**2))-acos((C1-sx)/sqrt(C1**2+C2**2))))/Wc
или
dt1=abs(acos((C1)/sqrt(C1**2+C2**2))-acos((C1-sx)/sqrt(C1**2+C2**2)))/Wc

где sx -- расстояние до края ячейки

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

thunar ★★★★★
() автор топика
Ответ на: комментарий от thunar

NGP - нулевой порядок точности. Попробуй лучше CIC - первый порядок точности. Насколько помню из-за NGP ошибки бысто накапливаются. Что с heating and deflection times? И насчёт дискретизации времени ЕМНИП для устойчивости leapfroga нужно, чтобы omega_p*delta_t<2.

Rakot ★★
()
Ответ на: комментарий от thunar

А эти самые колебания быть не должны, с точки зрения физики? Проверь. А то мало ли что...

TheKnight ★★★
()
Ответ на: комментарий от thunar

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

dmfd
()
Ответ на: комментарий от Rakot

лучше CIC - первый порядок точности

учитывая, что leap-frog имеет второй порядок --- в идеале надо и формфактор частиц брать соответствующим --- например треугольник. Это дает аппроксимацию полей по трем точкам.

нужно, чтобы omega_p*delta_t<2

это в теории. а на практике omega_p*delta_t<=0.2

anonymous
()
Ответ на: комментарий от anonymous

а вот как здесь http://www.dissercat.com/content/modelirovanie-dinamiki-chastits-v-darvinskoi... (http://ompldr.org/vZW43ZA/МОДЕЛИРОВАНИЕ ДИНАМИКИ ЧАСТИЦ В ДАРВИНСКОЙ ПЛАЗМЕ-d...) утверждают, что смогли делать расчёт при \omega_p\cdot\delta t > 2, но пока я об эти неявные схема только голову сломал.

thunar ★★★★★
() автор топика
Ответ на: комментарий от thunar

ну так вот если решаю уравнения движения в виде...

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

можно записать уравнения движения в виде...

Вопрос на засыпку: что будет, если в какой-то ячейке магнитное поле нулевое?

anonymous
()

В чём дело и почему нельзя применять полуаналитическое описание?

Главное --- потому что «полуаналитическое» описание предполагает постоянное внешнее электрическое поле. А Вы его на каждом шаге пересчитываете (следуя методу PIC).

Остальное детали.

Короче, делайте PIC высокого порядка, много частиц в ячейке, маленький шаг --- и будет Вам счастье. VLev

anonymous
()
Ответ на: комментарий от anonymous

потому что «полуаналитическое» описание предполагает постоянное внешнее электрическое поле

Я с самого начала об этом говорил, ня! ^____~

dmfd
()
Ответ на: комментарий от anonymous

схема Бориса

она и есть

если в какой-то ячейке магнитное поле нулевое?

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

thunar ★★★★★
() автор топика
Ответ на: комментарий от thunar

а вот как здесь

Знаю я этого Борадачёва. И модель Дарвина тоже знаю.

утверждают, что смогли делать расчёт при \omega_p\cdot\delta t > 2

Можно. Но чудес не бывает. Большой шаг -> нет плазменных колебаний. Тогда уж делайте как все --- перенормируйте диэлектрическую постоянную вакуума. VLev

anonymous
()
Ответ на: комментарий от dmfd

как сохранение полной энергии можно ввести в численную схему с многими частицами

Строго говоря, можно. Причем, разными способами. Я даже знаю один правильный. Но он редко используется, ибо не сохраняет полный импульс. Кстати, сам по себе метод PIC --- следствие законов сохранения. Но уравнений там хватает только на два из трех [физических] законов. обычно выбирают заряд и импульс.

anonymous
()
Ответ на: комментарий от thunar

появляется несколько электронов, способных перескочить 2 узла одновременно

Могу навскидку предложить целых два «метода», позволяющих «продолжить счёт» без особых проблем: 1. Уменьшать шаг по времени по мере появления быстрых частиц 2. «Разрешить» частицам перескакивать любое количество ячеек за шаг

Хуже от этого не станет --- точно. :)

anonymous
()
Ответ на: комментарий от anonymous

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

Большой шаг -> нет плазменных колебаний.

А есть ли формализм, позволяющий рассматривать плазму на временах t>>1/wp, считая её за это время стационарной без плазменных колебаний, но при этом сохранив кинетическое описание, так как оно необходимо для расчёта столкновений? на ум приходит только решать последовательными приближениями. Запомнить началные фазовые координаты, решить ур-я движения за время t, раздав плотности по узлам, решить пуассона. Потом снова прогнать решения со старыми фазовыми координатами, но пр новых полях. И делать итерации, пока не сойдётся.

thunar ★★★★★
() автор топика
Ответ на: комментарий от anonymous

Я даже знаю один правильный

Можно ключевое слово? Правда интересно. Всё моё теоретическое мировоззрение, основанное на скобках Пуассона и интегралах движения, против этого протестует.

dmfd
()
Ответ на: комментарий от thunar

в работе Батищева

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

А есть ли формализм...

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

стационарной без плазменных колебаний

это как раз из области чудес. Кстати, «плазменные колебания» --- прямое следствие уравнения Пуассона

при этом сохранив кинетическое описание

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

anonymous
()
Ответ на: комментарий от dmfd

Можно ключевое слово?

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

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

anonymous
()
Ответ на: комментарий от thunar

а так можно?

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

anonymous
()
Ответ на: комментарий от thunar

вот от этой работе речь

угу. Я просто ее читал 11 лет назад, соответственно. К тому же, про приложения тогда больше на Токамак обращал внимание. Кстати, про «перенормировку» в приложении В написано. Для частиц то же можно сделать, и даже проще IMHO. Только непонятно к чему это приведет в дальнейшем.

anonymous
()

это я вчера от anonymous-а писал

VLev
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.