Диссертация (1104059), страница 10
Текст из файла (страница 10)
Такжебыл описан метод построения решения с использованием частныхпроизводных по пространственным координатам от начальных данных.Следует отметить, что не для любых начальных распределений плотностизаряда удается получить аналитическое решение. Поэтому рассмотримпостроение численных методов для решения DV-задачи.2.1 Порядок аппроксимации решенияНа основании алгоритма, использованного для получения выражений((31)-(34), §1 глава 2) для коэффициентов разложений ((21), §1 глава 2),можно построить численный алгоритм необходимого порядка точности.Рассмотрим аппроксимацию из трех первых членов ряда ((21), §1 глава2).
Для каждого шага по времени получаем аппроксимацию решения в виде:t2D p, t D p,0 t Dt p,0 Dtt p,0 ,22tv p, t v p,0 t vt p,0 vtt p,0 .2(1)Коэффициенты при первой степени по времени t полностью выражены черезначальные условия по формулам ((31), (32) §1 глава 2), а коэффициенты((33), (34) §1 глава 2) при второй степени времени t 2 выражаются черезкоэффициенты при первой степени по времени. Поэтому выразимкоэффициенты ((33), (34) §1 глава 2) только через начальные условия,получим:96 2vp,0 V0 , D0 p V0 divD0 p rot V0 , D0 p .2 t000(2)Теперь получим соотношение для коэффициента Dtt p,0 2Dp,0 D0 p divD0 p V0 divD0 p ,V0 2 t0(3) rot V0 , rot V0 , D0 p .Длячисленного расчета можно использовать формулы(1) скоэффициентами, посчитанными по формулам ((31)-(34) §1 глава 2) или((31), (32) §1 глава 2) и (2), (3).
В качестве примера приведем расчет задачи(14) на двух схемах (1) с первым и вторым порядком аппроксимации повремени.2.2 Разностная схема с первым порядком аппроксимацииРанее было рассмотрено аналитическое решение задачи для однороднозаряженного шара. Опишем разностную схему для решения такой задачи.Пусть есть трехмерная область, в которой предполагается решать задачу. Дляопределенности в качестве геометрической формы области возьмемпараллелепипед с длинами сторон Lxs , s 1,2,3 .
В этой области зададимпрямоугольную сеткуi, j, k 1 i N x1 1, 1 j N x2 1, 1 k N x3 1 ,97(4)с шагом hxs LxsN xs, где N xs − число разбиений отрезка Lxs . По временизададим сетку 0 n NT 1 с шагом T, где T − промежуток времени, вNTкотором решается задача. Система уравнений ((18) §1 глава 2) в разностномвиде для первого порядка аппроксимации по времени примет вид: x1 , x2 , x3 R 3 ,Di, js,k ,n1 Di, js,k ,n i , j ,k ,nVi , j ,sk,n ,xxxs 1,2,3 xxxxVi , j ,sk,n1 Vi , j ,sk,n Di, js,k ,n dVi , j ,sk,n , 0Di1,1 j ,k ,n Di1,1 j ,k ,nxi , j ,k ,n x2hx1Di, j21,k ,n Di, j21,k ,nxx2hx2Di, j3,k 1,n Di, j3,k 1,nx(5)x2hx3,dVi , j ,sk,n Vi , j 1,k ,n w1i, js,k ,n Vi , j ,2k,n w2i, js,k ,n Vi , j ,3k,n w3i, js,k ,n ,xxгде выраженияxxxxw1i , js,k ,n , w2i, js,k ,n , w3i, js,k ,nxxxxявляются производными отскорости и определяются в соответствии с формулами:Vi 1,s j ,k ,n Vi 1,s j ,k ,nx xs w1i , j ,k ,n x2hx1w3i , js,k ,n xVi , j s 1, k ,n Vi , j s 1, k ,nx xs , w2i , j ,k ,n xs x2hx2(6) xs Vi , j ,k 1,n Vi , j ,k 1,n2hx3,.Поиск решения происходит следующим образом.
Сначала задается начальноеDi, js,k ,0xраспределениеi, j, k 0 i N x ,1Vi , j ,sk,0 ,xигде0 j N x2 , 0 k N x3 . Далее по формуле (5) сиспользованием формул (6) определяем Di, js,k ,1 и Vi , j ,sk,1 в узлах сетки (4).xxИспользуя граничные условия на поверхности S или условия симметрии98задачи, находим недостающие значения Di,xjs,k ,1 и Vi , xj ,sk,1 в узлах сеткиi, j, k i 0, N x ,1j 0, N x2 , k 0, N x3 . Аналогичным образом находимзначение сеточных функций Di, js,k ,n и Vi , j ,sk,n на следующих слоях при n 1 .xxПример 1Указанную задачу решим двумя способами.
Первый способ - путемиспользования разностной схемы (5), а второй методом «частица на частицу»PP (Particle-to-Particle). Оба результата сравним с известным аналитическимрешением ((44) §1 глава 2), полученным ранее.Итак, в силу симметрии задачи, воспользуемся сферической системойкоординат. В качестве начальных условий возьмем следующие величины:V0 r , D0 r er0 03r , r 0, R0 ,Q4, Q N p q , V R03 ,V3(7)R0 2m , N R 200, T 0.1sec,N p 2000, NT 200.Здесь использованы обозначения.
0 − начальная плотность заряда в шаре;Q − суммарный заряд шара; V − начальный объем шара; N p − числокрупных частиц для моделирования методом частица на частицу (PP); q −заряд одной крупной частицы; R0 − начальный радиус шара; N R − числоузлов разностной сетки вдоль радиуса; T − временной интервал, в течениекоторого происходит эволюция системы; NT − число шагов по времени.Результаты численных расчетов приведены на рис.1а и 1б. На рис.1апоказано распределение плотности заряда r вдоль радиуса. Синимцветом, показано распределение плотности, соответствующее задаче ((18) §199глава 2).
Красным цветом показана гистограмма, описывающая метод PP.Графики, приведенные на рисунке, соответствуют начальному и конечномумоменту времени. Приведено теоретическое решение (пунктирная линия).Рис.1б аналогичен рис.1а с той разницей, что на нем приведены графики r r 2 . Рис.1б дает распределение линейной плотности вдоль радиуса, и нанем более наглядно видно сохранение площади под кривой, котораясоответствует суммарному заряду.Из графиков на рис.1а,б видно, что разностная схема (5) дает хорошеесовпадение с теоретическим результатом. Для иллюстрации метода «частицана частицу» (PP) используется гистограмма. В начале координат при радиусеравном нулю аппроксимация функции плотности гистограммой не оченьуспешна, на рисунке видна некая осцилляция функции плотности.
Связаноэто с тем, что для нахождения плотности необходимо разделить друг надруга две маленькие величины: заряд на объем. Их малость вызвана тем, чтопри движении к началу координат, с уменьшением радиуса объемсферического слоя или шара уменьшается.
Следовательно, и количествозаряда, содержащееся в таком объеме тоже должно уменьшаться, так какплотность постоянна. В результате получается большая погрешность вгистограмме функции плотности при численном генерировании координатчастиц для метода PP.На рис.2 показана эволюция функции плотности заряда. Сплошнойлинией показан теоретический график, а пунктиром – расчетный. Здесь стоитотметить, что исходя из решения ((44) §1 глава 1) на каждом шаге по временивнутри шара плотность не зависит от радиуса и остается постояннойвеличиной, то есть зависит только от времени. Поэтому на рис.2 нетзависимости от координаты, то есть данное значение плотности можносопоставить любой точке внутри шара. Как видно из графика на рис.2 имеетместо хорошее совпадение между теорией и решением, полученным по схеме(5), которая имеет всего лишь первый порядок аппроксимации по времени.100Рис. 1а,б.
Начальное и конечное распределение объемной плотности частицдля модели шара с постоянной плотностью.Рис. 2 Эволюция объемной плотности частиц от времени для модели шара спостоянной плотностью.На рис.3 показано конфигурационное пространство для метода«частица на частицу» (PP). Слева дано распределение частиц в начальный101момент времени, а справа показано конечное положение частиц. Из рисункавидно объемное расширение шараРис. 3. Конфигурационное пространство начального и конечногораспределения частиц для модели шара с постоянной плотностью.Пример 2Рассмотрим случай, когда распределение заряда неоднородно вдольрадиуса.
Пусть функция плотности заряда имеет вид:n r 12 re ln r 222, r qN p n 2r ,2 r 2(8)где n r − нормальное логарифмическое распределение; N p − числочастиц. В качестве начальных условий возьмем следующие значения:N p 1000, 0, 0.2m ,reV0 r , D0 r r2 s 2 s ds,r 0r 0,1 m , N R 200, T 0.01sec, NT 400.102(9)На рис. 4а и 4б приведено начальное и конечное распределениефункции плотности вдоль радиуса. Сплошной линией показана функцияплотности, полученная по разностной схеме (5). В виде гистограммыприведена плотность частиц для расчета по методу «частица на частицу».
Нарис. 4а показана функция r , а на рис.4б функция r r 2 .Рис. 4. . Начальное и конечное распределение объемной плотности частицдля модели шара с нормальным логарифмическим распределениемплотности заряда.Из сравнения распределений на рис.4 видно, что новая постановка иметод PP дают похожий характер эволюции функции плотности.На рис.
5 показана эволюция функции плотности заряда через равныепромежутки времени для постановки задачи ((18) §1 глава 2), здесь видно,как происходит размывание сферического слоя заряда.103Рис. 5. Эволюция объемной плотности частиц от времени для модели шара снормальным логарифмическим распределением плотности заряда.Аналогичный расчет был проведен для трехмерного случая по схеме(5) на прямоугольной сетке:N x1 N x2 N x3 100, NT 20.(10)На рис.
6 показано пространственное распределение функцииплотности в начальный момент в медианной плоскости, а на рис. 7 показаноаналогичное распределение через момент времени T .104Рис. 6. Начальное распределение объемной плотности заряда в медианнойплоскости для модели шара с нормальным логарифмическим распределениемплотности зарядаРис. 7. Конечное распределение объемной плотности заряда в медианнойплоскости для модели шара с нормальным логарифмическим распределениемплотности заряда105Пример 3Рассмотрим еще один интересныйслучай, а именно примеравтомодельного решения, описанного выше. Для этого в примере номер 2добавим начальную постоянную скорость V0 . Исходя из общих соображений,изложенных выше, в этом случае должно получиться то же самое решение,что и в примере 2, только движущееся как целое со скоростью V0 .Как и в предыдущем случае, будем решать трехмерную задачу напрямоугольной сетке. В силу того, что задача теряет симметрию, величинаV0 , D будет отлична от нуля, и необходимо решать полную постановкузадачи ((14) §1 глава 2).Для численного решения сначала воспользуемся аппроксимациейпервого порядка по времени.















