2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Радиационное трение
Сообщение02.04.2023, 21:05 


04/03/23
12
Здравствуйте! Существует известное уравнение движения частицы, учитывающее радиационное трение.
$m\ddot{r(t)}=F_{rad}+F_{ext}(t)=\frac{q^2}{6\pi\varepsilon_{0}c^3}\dddot{r(t)}+F_{ext}(t)$
И если считать, что $F_{ext}=0$, то возникают так называемые "runaway solutions"
$\ddot{r}=\tau\dddot{r}\   \Rightarrow   r(t)=r_{0}+v_{0}t+a_{0}e^\frac{t}{\tau}$
А если, допустим, поместить частицу в постоянное электромагнитное поле, т. е $F_{ext}=q[E+v(t)\times B]$,
$m\ddot{r(t)}=\frac{q^2}{6\pi\varepsilon_{0}c^3}\dddot{r(t)}+q[E+v(t)\times B]$
то кажется, что движение частицы по спирали будет затухать.
Я хотел получить графики зависимостей: $v(t), a(t)$
Однако, численно решая соответствующие диффуры(опыта у меня в этом немного), я все равно получаю бесконечные runaway решения, а не ожидаемое затухание. Возможно, моя ошибка заключается в неверном подборе начальных условий. Прошу помощи, если понимаете какие начальные условия можно принять. Возможно,у меня есть ошибки в рассуждениях

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение03.04.2023, 09:44 


17/10/16
4757
AlexFuser
По моему, "runaway solution" - это просто обращенное во времени, нефизичное по сути, решение, т.к. сила радиационного трения тут взята с противоположным знаком (заряд не излучает и тормозится, а поглощает излучение и ускоряется). Такое же решение, как и обращение во времени решения для механической системы с трением. Следует рассматривать уравнение вида: $\ddot r=-\tau\dddot r$
Аналогично и для второго случая.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение03.04.2023, 21:29 
Заслуженный участник


03/01/09
1700
москва
Сила радиационного трения мала по сравнению с силой Лоренца, поэтому, мне кажется , применим метод последовательных приближений. В нулевом приближении считаем $F_{rad}=0$ и находим траекторию движения заряда $\vec  r_0(t)$. Затем в уравнении 1-ого приближения в выражение для $F_{rad}$ подставляем $ \dddot r_0(t)  $ и т.д.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение03.04.2023, 23:17 


30/01/23
17
Непростое это дело.
Академик пишет https://www.mathnet.ru/links/3f2bd92151 ... n11203.pdf

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение03.04.2023, 23:51 
Заслуженный участник


18/09/21
1756
ТС спрашивал только про вычислителный подход.
Соглашусь, с mihiv, что тут достаточно первой поправки (сама формула радиационного трения тоже является приближением).
Сначала считаем траекторию без этой силы. Потом считаем траекторию с этой силой определяемой по старой траектории.
Само собой, сравниваем, что расхождение между двумя траекториями мало.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение04.04.2023, 01:03 
Заслуженный участник
Аватара пользователя


23/07/08
10887
Crna Gora
Попробуйте, исходя из дифференциального уравнения и конечных условий $\mathbf r(t_0),\dot{\mathbf r}(t_0), \ddot{\mathbf r}(t_0)$, строить решение от момента $t_0$ назад во времени. Понятно, что это несколько другая задача, и на какие-то Ваши вопросы она не ответит; но на какие-то всё же ответит.
Вопрос, откуда взять конечное ускорение $\ddot{\mathbf r}(t_0)$, не так важен — хотя бы из условия $\dddot{\mathbf r}(t_0)=0$. Ошибка от неправильного задания ускорения будет сказываться только в окрестности $t_0$, в том смысле, что если эту окрестность отбросить, остальная часть решения будет одним из правильных движений частицы в присутствии силы радиационного трения.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение04.04.2023, 22:53 


05/12/21

138
AlexFuser насколько помнится, радиус кривизны (до релятвиской поправки) не зависит от скорости, но только от магнитного поля, заряда и массы частицы, т. е. заряженная частица будет двигаться по одной окружности, а с потерей энергии с разным периодом. Поэтому считать надо по углу поворота, а время просто суммировать, пока частица не достигнет минимальной скорости, предположим $1.0$ мм/сек :roll:

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение05.04.2023, 00:08 
Заслуженный участник
Аватара пользователя


23/07/08
10887
Crna Gora
Кстати, случай, когда поле $\mathbf E=0$, поле $\mathbf B$ однородно и частица движется в плоскости, перпендикулярной $\mathbf B$, решается аналитически.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение05.04.2023, 08:18 


17/10/16
4757
LLeonid3
Сила Лоренца пропорциональна скорости, а центростремительное ускорение при фиксированном радиусе - квадрату скорости. Так что с ростом скорости радиус окружности возрастает.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение05.04.2023, 11:13 


05/12/21

138
sergey zhukov да у меня заскок случился, думал о периоде, а написал радиус :oops:
Это период вращения заряженной частицы НЕ зависит от её скорости, а радиус "орбиты" как раз и меняется от скорости частицы (энергии), поэтому и числено интегрировать надо по углу, а не по времени, погрешность получится на порядки меньше.
Программку что-ли сделать :roll:

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение07.05.2023, 19:20 
Заслуженный участник


29/09/14
1239
Продолжение из этой темы. Рассмотрим частный пример движения частицы с зарядом $q$ и массой $m:$ пусть $\mathbf{E}=0,$ частица движется в плоскости, перпендикулярной вектору постоянного однородного магнитного поля $\mathbf{B}.$

Орт декартовой системы координат $\mathbf{e}_z$ выберем вдоль $\mathbf{B}.$ Тогда $\mathbf{B}=B\mathbf{e}_z,$ радиус-вектор частицы есть $\mathbf{r}(t) = x(t)\mathbf{e}_x+y(t)\mathbf{e}_y.$

Уравнение Ньютона (в гауссовой системе единиц) $$m \ddot{\mathbf{r}} = \dfrac{q}{c}[\dot{\mathbf{r}} \times \mathbf{B}]+\dfrac{2q^2}{3c^3} \dddot{\mathbf{r}}$$ сводится к системе двух уравнений для координат $x(t)$ и $y(t)$ частицы: $$\ddot{x}=\omega_c \dot{y} + \tau \dddot{x}$$ $$\ddot{y}=-\omega_c \dot{x} + \tau \dddot{y},$$ где $$\omega_c = \dfrac{qB}{mc}\,, \qquad \tau = \dfrac{2q^2}{3mc^3}\,.$$

При $\tau=0$ частица вращалась бы по окружности постоянного радиуса с "циклотронной частотой" $\omega_c.$ Четыре постоянных интегрирования (начальные данные): координаты центра орбиты $x_c$ и $y_c$, начальная фаза вращения $\alpha$, радиус орбиты $C.$ Полагая для удобства, что начало координат помещено в центр орбиты (так что $x_c=y_c=0),$ можно один раз проинтегрировать уравнения с равными нулю постоянными интегрирования; такие же уравнения получатся для переменных $x-x_c$ и $y-y_c$ с произвольными $x_c,\,y_c.$

Поступим так же при $\tau > 0.$ Вот получившиеся дифференциальные уравнения (ДУ) в записи с безразмерной переменной времени $t'=\omega_c t$ и с безразмерным неотрицательным параметром радиационного трения $w=\omega_c \tau$ $$\dot{x} = y +  w\ddot{x},$$ $$\dot{y} = -x + w\ddot{y}.$$
1. Сначала найдём точное решение, чтобы было с чем сравнивать приближённые ответы. "Методом комплексных амплитуд" ищем решение в виде $$x = \operatorname{Re}\left (A\,\exp(i\Omega t')\right ),$$ $$y = \operatorname{Re}\left (C\,\exp(i\Omega t')\right ),$$ где комплексные постоянные $A,\,C,\,\Omega$ определяются системой уравнений $$(\Omega^2w+i\Omega)A-C=0$$ $$A+(\Omega^2w+i\Omega)C=0.$$ Такая система имеет не равные нулю решения $A,C$ только тогда, когда её определитель равен нулю: $$(\Omega^2w+i\Omega)^2+1=0$$ Решение этого уравнения для $\Omega$ можно выбрать так, что: $$\Omega^2w+i\Omega=i,\qquad A=-iC,$$ $$\Omega = \dfrac{i}{2w}\left ( -1 +\sqrt{1-i4w} \right )$$ Обозначим: $G$ - мнимая часть, $R$ - вещественная часть частоты $\Omega,$ т.е. $$\Omega = R + iG$$ Положим $C=1$ - это означает, что начальную фазу $\alpha$ выбираем равной нулю, а расстояние от начала координат до частицы при $t=0$ выбираем за единицу длины. Тогда $A=-i.$

Таким образом, вот точные выражения для траектории частицы c $\alpha=0,\,C=1:$ $$x(t')=e^{-Gt'}\sin(Rt'),$$ $$y(t')=e^{-Gt'}\cos(Rt').$$ Для перехода к исходным единицам достаточно умножить $x$ и $y$ на произвольный вещественный коэффициет $C$ с размерностью длины и подставить $t'=\omega_c t,$ где $\omega_c=\frac{qB}{mc}.$ К аргументам синуса и косинуса можно добавить произвольное вещественное число $\alpha.$ Если к получившимся $x(t)$ и $y(t)$ добавить произвольные значения координат центра орбиты $x_c$ и $y_c,$ то в ДУ следует под координатами $x$ и $y$ понимать разности $x-x_c$ и $y-y_c.$

В решении $\Omega = R+iG$ уравнения $\Omega^2w+i\Omega=i$ величина $R$ имеет смысл частоты вращения, а $G$ - параметр "затухания спирали"; оба параметра безразмерные, измеряются в долях от $\omega_c.$ Из уравнения для $\Omega$ следует, что $w(R^2-G^2)-G=0$ и $2wRG+R=1.$ Видно, что в отсутствие трения, т.е. при $w=0$ получается, как и должно быть, $R=1,\,G=0,$ а при небольшом трении, т.е. при $0<w\ll 1$ получается немного меньшее единицы значение $R$ и маленькое положительное $G.$ Значит, с трением спираль будет не раскручиваться, а сжиматься, "затухать".



2. Теперь для упражнения попытаемся решать ту же систему ДУ приближённо.

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

Переходим к системе ДУ первого порядка c вспомогательными функциями $p(t')$ и $r(t')$: $$\dot{x}=p,\quad \dot{p}=(p-y)/w, \quad \dot{y}=r, \quad \dot{r}=(r+x)/w.$$ Шаг по безразмерному времени $t'$ обозначаю как $h.$ Приращения функций записываю в виде оборванного ряда по степеням $h$ с точностью до членов порядка $h^4$ включительно; например, для функции $p$ имеем приближённое равенство: $$p(t'+h)=p(t')+\dot{p}h+\frac{1}{2}\ddot{p}h^2+\frac{1}{6}\dddot{p}h^3+\frac{1}{24}\ddddot{p}h^4,$$ где производные относятся к точке $t'.$ Затем эти производные последовательно выражаем через сами функции в точке $t'$ с помощью наших ДУ; например, $\dot{p}=(p-y)/w,$ дифференцируя это ещё раз и снова пользуясь ДУ, получаем выражение для $\ddot p,$ и так далее.

Значения функций в точке $t'$ пишем с индексом, например, $n.$ В точке $t'+h$ они тогда будут с индексом $n+1.$ То есть: $p_{n+1}=p_n +...,$ и так далее. Так получается система рекуррентных соотношений, позволяющих в цикле по $n$ от нуля до какого-нибудь большого $N\gg 1$ вычислять интересующее нас множество значений координат $x_n,\,y_n,$ вместе со вспомогательными $p_n, \, r_n,$ имеющими смысл значений скорости $\dot{x}_n,\,\dot{y}_n.$

a) Посмотрим для начала, что получается при положительном шаге $h=0.01$ с параметром трения $w=0.2.$

Пусть начальные координаты частицы: $x_0=0,\, y_0=1.$ При этом из приведённых выше формул точного решения (c $\alpha=0,\, C=1)$ видно, что начальную скорость можно выбрать так: $p_0=\dot x(0) = R,$ $r_0=\dot y(0) = -G,$ где $R$ и $G$ определяются из решения уравнения для $\Omega=R+iG.$ Вычислить $\Omega$ поручаю Маткаду (но не знаю, насколько хорошо он с таким заданием справляется :). Вот маткадные фрагменты расчёта, рекуррентные формулы:

Изображение

Ниже на графиках красная линия - точное решение, пунктир - результат приближённого вычисления по указанным рекуррентным формулам. Пунктир это способ изображения графика, а не сами вычисленные точки. Приближённых точек так же много, как и точных ($N=2000),$ но если их все изобразить, то они полностью закрашивают собой предыдущий график. Результат при $h=0.01$ и $w=0.2:$

Изображение

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

б) Выполняем расчёт "назад во времени" - при отрицательном $h=-0.01,$ с прежним значением параметра трения $w=0.2.$

Начальные значения координат и скорости в этом примере выбираю те, которые получаются из точного решения для конечного момента времени $t'=20,$ т.е. для точки с номером $N=2000:$ $x_0=-4.0379830047\cdot 10^{-3},$ $y_0=3.3374190164\cdot 10^{-2},$ $p_0=3.1938495138\cdot 10^{-2},$ $r_0=-1.88\cdot 10^{-3}.$ После завершения цикла вычислений (с прежними рекуррентными формулами, по $n$ от нуля до $N=2000,$ но теперь с отрицательным шагом $h<0),$ все получившиеся координаты $x_n,\, y_n$ перенумеровываем в обратном порядке. И... ура, неустойчивость исчезла, результат соответствует точному решению:

Изображение

"Затухающая спираль" на плоскости $x,y$ в этом примере выглядит так:

Изображение

в) Расчёт "назад во времени", при $h=-0.01,$ с прежним значением параметра трения $w=0.2,$ но с произвольными начальными условиями:

Начальные условия теперь задаю как попало; например, пусть $x_0=y_0=p_0=r_0=0.1.$ По отношению к нумерации "вперёд во времени" эти условия конечные - они для точки $t'$ с номером $N=2000.$ После завершения такого же, как и раньше, цикла вычислений по $n$ от нуля до $N=2000,$ с отрицательным шагом $h=-0.01,$ получившиеся координаты $x_n,\, y_n$ перенумеровываем в обратном порядке. В результате, получившиеся $x_0,\,y_0$ относятся уже к точке $t'=0.$ По ним определяем начальный радиус $C=\sqrt{x_0^2+y_0^2}$ и начальную фазу $\alpha=\arctg \frac{x_0}{y_0}$ для сравнения с точным решением $$x(t')=Ce^{-Gt'}\sin(Rt'+\alpha),$$ $$y(t')=Ce^{-Gt'}\cos(Rt'+\alpha).$$ Оказывается, в этом случае результат приближённого вычисления тоже соответствует точному решению:

Изображение

Лишь на первых шагах расчёта "назад во времени" (т.е. вблизи $t'=20$ в этом примере) заметно обусловленное произвольными начальными условиями большое отличие приближённого результата от точного решения:

Изображение


г) Расчёт "назад во времени" при $h=-0.01$ с малым значением параметра трения: $w=0.01.$

Пусть максимальное $t'$ будет равно $1000,$ так что $N=1000/|h|=100000.$ Начальные условия (т.е. условия для точки $t'$ с номером $N)$ задаю произвольно: $x_0=y_0=p_0=r_0=0.00033.$ Сравнение с точным решением делается тем же способом, что и в предыдущем примере; только теперь $N$ в пятьдесят раз больше. Цель этого теста - убедиться, что даже на длинном куске траектории приближённый расчёт не сбивается, согласуется с точным решением. Оказывается, так и есть. Вот часть траектории в начале временного интервала:

Изображение

Последняя часть траектории:

Изображение

Вся "спираль" на плоскости $x,y$ в этом примере выглядит так:

Изображение

Маткад сумел даже кое-как анимировать примерно 80 оборотов по этой траектории (здесь ссылка на avi-файл, некоторые браузеры могут его показывать, но не знаю, все ли ).


О перенумеровке точек траектории в обратном порядке:

Есть варианты. Вычисленные по указанным в начале рекуррентным формулам (c $h<0)$ массивы $x_n, y_n$ можно в цикле по $j$ от нуля до $N$ переписать в новые массивы, например, $x1_j:=x_{N-j}, \, y1_j:=y_{N-j},$ а затем, чтобы вернуться к исходным обозначениям, в отдельном цикле записываем их на место исходных величин: $x_j:=x1_j$ и т.д. Более простой вариант: в рекуррентных формулах можно дать массивам новые обозначения, например, $x1_n,\,y1_n$ вместо $x_n, y_n$ и т.д.; тогда для перенумеровки в обратном порядке достаточно в цикле по $j$ от нуля до $N$ записать $x_j:=x1_{N-j},\,y_j:=y1_{N-j}.$ Наконец, самый простой вариант - вести рекуррентные вычисления сразу в обратном порядке, т.е. в цикле по $n$ от $N$ до $1$ с шагом $-1$ вычислять $x_{n-1}=\text{выражение с индексом n}$ и т.д.; тогда не надо будет ничего перенумеровывать.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение07.05.2023, 23:29 
Заслуженный участник
Аватара пользователя


23/07/08
10887
Crna Gora
Cos(x-pi/2)
Спасибо за Ваш труд. :-)

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение08.05.2023, 00:31 
Заслуженный участник


29/09/14
1239
svv
Огромное спасибо Вам!

Мне не за что (я лишь написал про то, что попытался понять на своём ученическом уровне из Ваших сообщений :-) :oops: ).

Все идеи - в Вашем потрясающем сообщении:
svv в сообщении #1588186 писал(а):
Попробуйте, исходя из дифференциального уравнения и конечных условий $\mathbf r(t_0),\dot{\mathbf r}(t_0), \ddot{\mathbf r}(t_0)$, строить решение от момента $t_0$ назад во времени. Понятно, что это несколько другая задача, и на какие-то Ваши вопросы она не ответит; но на какие-то всё же ответит.

Вопрос, откуда взять конечное ускорение $\ddot{\mathbf r}(t_0)$, не так важен — хотя бы из условия $\dddot{\mathbf r}(t_0)=0$. Ошибка от неправильного задания ускорения будет сказываться только в окрестности $t_0$, в том смысле, что если эту окрестность отбросить, остальная часть решения будет одним из правильных движений частицы в присутствии силы радиационного трения.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение08.05.2023, 18:58 


04/03/23
12
Cos(x-pi/2), я разобрался, спасибо! Я Вам очень благодарен!
svv, огромное спасибо!
Для меня это было очень полезно! Еще раз большое спасибо.

 Профиль  
                  
 
 Re: Радиационное трение
Сообщение08.05.2023, 19:42 
Заслуженный участник
Аватара пользователя


23/07/08
10887
Crna Gora
AlexFuser, Вам тоже спасибо за интересную задачу!

(Оффтоп)

Думал, не изложить ли своё точное решение с ложечкой дёгтя, но не уверен, что это нужно. :-)

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 21 ]  На страницу 1, 2  След.

Модераторы: photon, whiterussian, profrotter, Jnrty, Aer, Парджеттер, Eule_A, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group