close

Вход

Забыли?

вход по аккаунту

- Вестник Пермского университета. Математика

код для вставкиСкачать
ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА
Математика. Механика. Информатика
2014
Вып. 3(26)
УДК 537.2
Задача электростатики о взаимодействии
заряженных шаров на близких расстояниях
Е. Л. Тарунин
Пермский государственный национальный исследовательский университет
Россия, 614990, Пермь, ул. Букирева, 15
(342) 239-64-09
Известно, что сила взаимодействия заряженных проводников на близких расстояниях существенно отличается от закона Кулона для точечных зарядов, размещенных в центре шаров. Это отличие для проводящих сфер было подробно исследовано в работах [1–3], для тел
цилиндрической формы – в [5]. В работах [1–3] расчеты были выполнены с использованием
емкостных коэффициентов, вычисляемых с помощью бесконечных рядов [4]. В [5] вычисления были проведены на основе решения методом сеток задачи Дирихле для потенциала
электростатического поля. В работе [6] для вычисления сил взаимодействия использовано
решение уравнения Лапласа в бисферических координатах. В данном исследовании расчеты выполнены для проводящих сфер, как и в [1–3], но с использованием вычисляемого потенциала электростатического поля. Они показали как соответствие с результатами работ
[1–3], так и отличия.
Ключевые слова: электростатика; уравнение Лапласа; потенциал электростатического
поля; взаимодействие заряженных тел на близких расстояниях.
1. Постановка задачи
татов решения вблизи сфер использовались
сферические координаты. Для получения граничных значений для внешних радиусов сфер
была применена процедура интерполяции, неизбежно вносящая погрешность. Поэтому было решено вести счет во всей области, применяя уравнение Лапласа в цилиндрической системе координат. При этом возникают сложности счета в нестандартных узлах вблизи сфер.
И все же этот вариант расчета, обсуждаемый
ниже, признан наиболее оправданным.
Варианту алгоритма, результаты которого будут обсуждаться в этой статье, предшествовали два варианта. В итоге от них
пришлось отказаться. Все же кратко упомянем их.
В первом варианте счет выполнялся с
использованием уравнения Лапласа, записанного в цилиндрических и сферических координатах. Сферические координаты использовались для счета вблизи обеих сфер. Соответствующие области перекрывались, и для организации счета требовалось согласование
значений на их границах. От этого варианта
пришлось отказаться, так как было обнаружены затруднения в выборе областей при малых
расстояниях между шарами, особенно для
сфер с различными радиусами.
Во втором варианте основной счет выполнялся с использованием уравнения Лапласа, записанного в цилиндрической системе
координат. Для облегчения обработки резуль-
Геометрия расчетной области изображена на рис. 1. Расчеты выполнялись с учетом
осевой симметрии. Поэтому на рисунке изображена лишь половина области.
Рис. 1. Геометрия расчетной области
© Тарунин Е. Л., 2014
16
Задача электростатики о взаимодействие заряженных шаров …
Параметрами задачи являются расстояние между центрами сфер L=x2 - x1, радиусы
сфер R1, R2 и заряды на сферах Q1, Q2. Расчеты
выполнялись в безразмерных переменных, в
качестве единицы расстояния был выбран радиус первого шара R1 =1 (полагалось, что
R2  R1 ). Постановка задачи предполагает,
что взаимодействующие заряды изолированы.
Положение внешней границы, на которой задавалось нулевое значение потенциала, определялось параметром метода  >10‫׃‬
x1= AO1=  +R2, O2C=  +R2,
АB=CD=  +R2.
(1.1)
Уравнение Лапласа для потенциала в
цилиндрических координатах x, r с учетом
симметрии имело вид
1    2
 (r )  2  0 .
r r r
x
Постановка задачи позволяла также определить электроемкость шаров. Согласно
электростатике [9] в системе проводников заряды на проводниках связаны с потенциалами
на них соотношением
qi   ci , j j .
(1.3)
j
Здесь
а
ci , j
ci ,i
– коэффициенты емкости i-го тела,
– коэффициент электростатической ин-
дукции между телами с индексами i,j (этот
коэффициент определяет величину заряда qi,
когда потенциал j-го тела равен  j , а все остальные проводники заземлены). Для коэффициентов в формуле (1.3) выполняются соотношения
ci ,i  0, ci , j  c j ,i  0 при i  j
(1.2)
(1.4)
В случае двух проводников уравнения
(1.3) имеют вид
В стандартных узлах уравнение аппроксимировалось
на
квадратной
сетке
r  x  h . Соответствующая система уравнений решалась итерационным методом. Особая аппроксимация уравнения (1.2) использовалась на оси и в узлах с нестандартными шаблонами вблизи шаров. Описание этих аппроксимаций дано при описании метода.
При заданных значениях потенциалов
на сферах Vs1 и Vs2 из решения задачи Дирихле находилось распределение потенциала.
При равных значениях потенциалов на сферах
Vs2=Vs1 решалась задача о взаимодействии
одноименных зарядов. Задание Vs2=-Vs1
приводило к задаче о взаимодействии разноименных зарядов.
По вычисленным значениям потенциалов находились различные характеристики
решения: сила электростатического взаимодействия F, зависимость радиальной компоненты поля на сферах Er от полярного угла
 , отношение максимальной радиальной
компоненты поля на сфере к минимальной
компоненте kE1 и другие.
Основной характеристикой решения являлось отношение вычисленной силы к силе
взаимодействия
по
закону
Кулона
kf  F0 / F . Величина kf показывает отклонение силы взаимодействия от закона Кулона для точечных зарядов (при kf  1 сила
q1  c1,1  1  c1,2   2 ,
q2  c2,1  1  c2,2   2 .
(1.5)
Коэффициенты в этих формулах зависят
от радиусов сфер и от расстояния между сферами. При равенстве радиусов (именно этот
случай рассматривается в статье) c1,1  c2,2 .
При задании равных потенциалов  2  1
реализуется вариант отталкивания одноименных зарядов, а при  2  1 – вариант притяжения зарядов разного знака. В общем случае
произвольных значений потенциалов реализуется вариант с любым отношением зарядов
  q1 / q2 
c1.1  c1,2 ( 2 / 1 )
. (1.6)
c1,2  c2,2 ( 2 / 1 )
Решение сформулированной задачи с
потенциалами 1  1,  2  0 и вычисление
зарядов позволяют из формул (1.5) определить коэффициенты
c1,1  q1 , c1,2  q2 .
(1.7)
Вычислительные эксперименты позволили выяснить, что при стремлении расстояния между сферами к бесконечности значения
q2 , c1,2 стремятся (как и положено) к нулю, а
взаимодействия меньше F0, а при kf  1
больше).
17
Е. Л. Тарунин
значения
q1 , c1,1
Для сокращения объема вычислений
значения компонент массивов В1 и В2 вычислялись один раз перед началом итераций.
Система уравнений для потенциала
электростатического поля решалась методом
последовательной верхней релаксации [7, 8].
Для узлов со стандартным шаблоном параметр
верхней релаксации вычислялся по формуле,
подобранной экспериментальным путем,
в выбранных единицах из-
мерения – к 4 (при L=6, например, отличие
q1 , c1,1 от 4 составляет 8%). Вычисленные
коэффициенты позволяют найти отношение
зарядов
  q / q 
Здесь
q
c1,1  c1.2
 1.
c1,1  c1,2
(1.8)
соответствует случаю одно-
именных зарядов ( 1   2  1 ), а
q

  2 /(1  sin(0.9 / NMR )) . (2.4)
– разно-
В этой формуле NMR – максимальное
значение индекса k по радиусу.
Для узлов с нестандартным шаблоном
вблизи сфер параметр релаксации равнялся 1
(или был чуть больше 1).
Первое слагаемое в уравнении Лапласа
(1.2) на оси имеет устранимую неопределенность типа ноль деленный на ноль. Применение к этому слагаемому правила Лопиталя
позволяет выяснить, что
именных ( 1   2  1 ). При удалении зарядов   1 , а c1,2  0 .
Основными параметрами сетки служили число интервалов на первом радиусе RN и
параметр  , определяющий внешний размер
области (см. (1.1)). Максимальное значение
индекса k равнялось NMR, а максимальное
значение индекса m равнялось NM. Для хранения всех элементов массива для потенциала
V[k,m] требовалось
1  
 2
lim r 0 (  (r ))  2( 2 ) (2.5)
r r r
r
n=RN2(L+2(R2+  ))(R2+  )
ячеек оперативной памяти. При типичных
значениях параметров задачи и метода (L=3,
RN=30,  =24, R2=1) число элементов массива для потенциала электростатического поля n
более миллиона.
и использовать с учетом симметрии аппроксимацию
4(V1,m  V0,m )
 2
.
2( 2 ) 
r
h2
2. Метод решения
(2.6)
Полная аппроксимация уравнения Лапласа на оси с учетом (2.6) приобретает вид
Аппроксимация уравнения Лапласа
на
квадратной
сетке
rk  k  h, xm  m  h в обозначениях школы
А.А. Самарского [7] имела вид
V0,m  (4V1,m  V0,m1  V0,m1) ) / 6 . (2.7)
(1.2)
Для сокращения общего числа итераций
значения потенциала на оси V0,m находились
1
((rk  h/2) Vr (rk  h/2) Vr ) Vx,x  0. (2.1)
rk
методом скалярной прогонки в предположении, что значения V1,m известны. Коэффициенты скалярной прогонки Pm , Qm вычислялись по рекуррентным формулам
Замена   V сделана из методических соображений для отличия разностного
решения от точного решения уравнения (1.2).
После введения обозначений
Pm 
B1[k] 0.25(10.5/ k), B2[k]  0.25(10.5/ k) (2.2)
1
, Qm  Pm (4V1,m  Qm1) (2.8)
6  Pm1
Прогонки осуществлялись для трех участков на оси (перед первой сферой, между
сферами и после второй сферы) в предположении, что значения потенциала на концах
участков заданы.
уравнение (2.1) записывалось в виде, удобном
для применения итерационного метода
V[k,m] (B1[k]V[k 1,m] B2[k]V[k 1,m]
(2.3)
V[k,m1]V[k,m1])/4.
18
Задача электростатики о взаимодействие заряженных шаров …
дартным шаблоном, наиболее удаленным от
поверхности соответствующей сферы. Первые
два вектора относятся к первой сфере, а следующие два – ко второй.
Итерации в стандартных узлах выполнялись с использованием подпрограммы с
именем "CalcV" в общем случае для 5 областей. Опишем границы этих областей, напомнив, что RN – число интервалов на первом
радиусе R1, RN2 – число интервалов на радиусе R2, NMR – номер граничного узла по вертикали, а NM – максимальный номер узла по
горизонтали.
Первая область содержит узлы перед
первой сферой, индекс k меняется в ней от 1
до RN, а индекс m – от 1 до mk1[k]-1.
Вторая область содержит узлы между
сферами, индекс k меняется в ней от 1 до RN,
а индекс m – от mk2[k]+1 до mk3[k]-1.
Третья область содержит узлы правее
второй сферы, индекс k меняется в ней от 1 до
RN2, а индекс m – от mk4[k]+1 до NM-1.
Четвертая область существует лишь при
выполнении неравенства R2>R1, индекс k меняется в ней от RN+1 до RN2, а индекс m – от
1 до mk3[k]-1.
Пятая область является наибольшей по
размеру, индекс k меняется в ней от RN2+1 до
NMR-1, а индекс m – от 1 до NM-1.
Компоненты векторов mk1[k], mk2[k],
mk3[k], mk4[k] вычислялись перед началом
основного счета с помощью алгоритма, названного условно – "лесенка". В качестве
примера опишем алгоритм вычисления компонент вектора mk2[k], предполагая, что нам
известно значение mk2[k-1] (при k=1, например, mk2[1]=M1+RN, где M1 – индекс m для
центра первой сферы). Алгоритм проще объяснять, используя вместо m значение dm[k],
равное отклонению индекса от оси сферы. В
цикле по k, начиная с k=1, проверяются условия нестандартности узла. Величины плеч
узлов пространственной сетки определятся в
масштабе длины стандартного узла h. Полагается, что найдено значение dm[k], если одно
из плеч узла менее 1, а для узла правее этого
узла все плечи равны 1.
Для реализации алгоритма требуется
вычислять расстояния от узла до поверхности
сферы по горизонтали w и вертикали s. Формулы для вычисления этих значений указаны
выше (2.10). При обработке результатов требуется вычислять радиальную компоненту
поля Er. Для нестандартного узла вблизи пер-
Рис. 2. Обозначения для узла с нестандартным шаблоном
Выпишем конечно-разностную аппроксимацию уравнения Лапласа для нестандартного шаблона, изображенного на рис. 2.
Для обозначения узлов использована
географическая символика с заглавными буквами N,E,S,W (N – северный узел, W – западный, E – восточный, S – южный; соответствующие маленькие буквы обозначают длину
плеч в единицах стандартного шага h).
Для узлов, расположенных слева от
сфер, w=1, а для узлов расположенных справа, e=1. В использованных обозначениях аппроксимация уравнения Лапласа имела вид
1
V V
V V
((rk h/2) N P (rk sh/2) P S )/(0.5(hsh))
rk
h
sh
(2.9)
1
V V V V
( E P  P W)0.
0.5(ehwh
 ) eh wh

Для использования итерационного метода это уравнение преобразовывалось так, чтобы можно было вычислять значение VP через
четырех соседей. Длины плеч для нестандартного узла с индексами k , dm вычислялись по
формулам (dm – номер индекса по горизонтали, отсчитываемый от центра сферы)
s  k  RN 2  dm 2 ,
w  dm  RN 2  k 2 .
(2.10)
Формулы приведены для узлов, расположенных справа от первой сферы. Подобный
вид имеют формулы для второй сферы и для
узлов, расположенных слева от сфер.
Наличие нестандартных узлов отражается и на организации циклов по узлам со
стандартным шаблоном. В программе вычислялись компоненты четырех векторов mk1[k],
mk2[k], mk3[k], mk4[k], которые определяли
значение индекса по горизонтали m с нестан-
19
Е. Л. Тарунин
вого шара радиальная компонента поля вычисляется по формулам
3. В этом блоке задаются нулевые значения
потенциала для узлов вне сфер и заданные значения потенциала Vs1, Vs2 внутри сфер и на их
поверхностях. Кроме того, задаются нулевые
значения для счетчика полных итераций
mfull=0. Полное число внешних итераций (циклов) обычно всегда было менее 200.
4. После сброса в ноль счетчика внутренних итераций m0 и значения невязки nev выполнялись итерации по внутренним узлам
расчетной области. Вначале по одному разу
вычислялись значения потенциала для узлов с
нестандартными шаблонами в следующем
порядке: 1) близнецы; 2) узлы при k=RN для
первого шара и k=RN2 для второго шара; 3)
счет в граничных нестандартных узлах для
обоих шаров с использованием обращения к
подпрограммам CaseLeft, CaseRight; 4) скалярные прогонки для трех участков на оси
r  0 . Далее осуществлялись итерации по
всем узлам со стандартным шаблоном. Эта
часть программы была наиболее трудоемкой.
В ней использовалось обращение к подпрограмме без параметров CalcV. Приведем текст
этой подпрограммы на Паскале АВС.
Procedure CalcV; Begin
dr  h k2  dm2  R1, Er  (V[k, m] Vs1)/ dr. (2.11)
Опишем используемую классификацию
нестандартных узлов. Назовем узлы с индексами mk1[k], mk2[k], mk3[k], mk4[k] для заданных значений k "граничными" нестандартными узлами (присвоение этим узлам названия граничных связано с тем, что рядом с ними слева или справа находится узел со стандартным шаблоном).
Анализ показывает, что для некоторых
значений k близких к k=RN есть нестандартные узлы, которые расположены ближе к
сфере. Такие узлы будем называть "близнецами". Число близнецов при фиксированном
k<RN может быть оценено по формуле, в которой использована процедура вычисления
целой части от вещественного числа:
max dm  trunc(2 RN 1  2  RN 1) .(2.12)
При RN=20 max dm =2, а при RN=100
max dm =5. Напомним, что число близнецов
указывается для узлов, расположенных с одной стороны сферы (общее число близнецов
вдвое больше).
Наибольшее число близнецов при
k=RN. Близнецы при k=RN выделены в отдельную группу. Максимальное число близнецов при k=RN равно
s1: ( B1[k ]  V [ k  1, m]  B 2[k ]  V [k  1, m] 
V [k , m  1]  V [k , m  1]) / 4  V [k , m];
V [k , m] : V [k , m]  OM  s1;
if abs(s1) >nev then nev:=abs(s1); end.
p1  trunc( 2  RN  1  0.000001). (2.13)
Заметим, что величина невязки nev вычислялась также в узлах с нестандартными
шаблонами. При завершении 4-го блока счетчик внешних итераций mfull увеличивался на 1
и выполнялось сравнение невязки с eps0=10-13.
Возврат в начало 4-го блока происходил при
одном
из
двух
условий
–
m0  NM , nev>eps0. Первое условие экономило время счета, сокращая число итераций на первых внешних циклах. Без первого
условия общее число внутренних итераций
Smo возрастало примерно на 40%. Обычно до
mfull  40 число итераций m0=NM, а затем
убывало до 1. Так, например, при
Отсюда следует, что при RN=10 p1=4, а при
RN=20 p1=6.
Перейдем к описанию основных блоков
программы.
1. В первом блоке описываются используемые переменные, подпрограммы и задаются значения параметров задачи и метода.
2. Во втором блоке вычисляются параметры нестандартных узлов: индексы граничных нестандартных узлов dmk1[k], dmk2[k],
показывающих удаление узла от оси сфер;
длины плеч нестандартных узлов hwk1[k],
hsk1[k] для первого шара и hwk2[k], hsk2[k]
для второго; p1, p2 – число нестандартных
узлов при k=RN и k=RN2 соответственно;
kd1[k[, kd2[k] – число близнецов. Заметим,
что часть массивов, упомянутых в этом блоке,
относятся лишь к правой стороне сфер. Для
левых половин соответствующие значения
вычисляются из соображений симметрии.
L  3,   10, RN  45, eps  1013
зависимость числа внутренних итераций mo
от числа внешних итераций mfull указана в
таблице:
mfull
mo
20
 30
1175
40
930
50
279
60
57
70
3
80
2
 90
1
Задача электростатики о взаимодействие заряженных шаров …
Обычно значения mo выводились с интервалом для внешних итераций, равном 10.
Для интегрального анализа сходимости на
печать выводились номера полных итераций
m1, m2 (m1 соответствует номеру, при котором впервые mo  NM , m2 соответствует
номеру, после которого все последующие
mo=1). Так, например, при значениях L=3,
 =30, RN=24 полное число итераций
Smo=73 424, а номера m1, m2 равны соответственно 40 и 102. Полное число внутренних
итераций Smo возрастало при увеличении
числа интервалов на первом радиусе RN
сильнее линейной зависимости (Smo  RN 1.2 ).
Для дополнительного контроля сходимости итерационного процесса через заданное
значение числа полных итераций на экран
компьютера выводилось значение суммы SSE:
CaseLeft, CaseRight; параметр верхней релаксации для узлов около сфер был близок к 1.
По вычисленным значениям потенциала
V[k,m] находились следующие характеристики решения: радиальные компоненты поля на
сферах En1[j], En2[j], заряды на полусферах
QL1, QR1, QL2, QR2, полные заряды на сферах Q1=QL1+QR1, Q2=QL2+QR2, силы f11,
f12,f21,f22, F1=f11+f12, F2=f21+f22. Значения
En1[j], En2[j] позволяли вычислить как заряды на сферах, так и компоненты сил. Интегральными характеристиками зависимости
En1[j] служили три величины – kE1 (перепад
напряженности на сфере), SEn: (среднее значение) и GladEn (характеристика гладкости):
kE1( j )  max( En1[ j ])/ min( En1[ j ]),
j
SEn 
2 JN 1
SSE (mfull ) 
 MEn1[ j ].
(2.14)
j
1
 En1[ j ].
2  JN1 j
(3.1)
Характеристика гладкости была введена
для того, чтобы сравнивать алгоритмы сглаживания. В первых расчетах эффекты негладкой зависимости напряженности от полярного
угла определялись по виду графического изображения функции En1(  ). Позднее была введена числовая характеристика GladEn. Для ее
вычисления в цикле по индексу j от 2 до
(2JN1-1) суммировались абсолютные значения выхода En1[j] за пределы интервала от
En1[j-1] до En1[j+1]. Затем в процентах вычислялось окончательное значение
j 0
Эта величина быстро стабилизировалась. Приведем в качестве примера изменения
этой величины при L=3.5,  =30, RN=24 через
2 шага. Относительная разница (SSE(4)SSE(2))/SSE(2) была равна 8.98%, а последующие разницы составляли 1.31 и 0.74%.
После 50 внешних итераций не изменялись
первые 8 значащих цифр SSE.
5. В этом блоке выполняется обработка
результатов решения. Для вычисления радиальных компонент поля En1[j], En2[j] вначале
вычисляются углы для узлов, ближайших к
поверхности сфер Mb1[j], Mb2[j], и соответствующие площади Mds1[j], Mds2[j] . По значениям радиальных компонент поля находятся отношения максимальной компоненты поля к минимальной kE1, kE2. Затем вычисляются заряды на половинках сфер QR1, QR2,
QL1, QL2 и силы взаимодействия f11, f12, f21,
f22. В конце обработки результатов расчета на
экране строится зависимость радиальной
компоненты поля от полярного угла. Завершается программа построением изолиний.
GladEn : 100GladEn / SEn /(2 JN1  1). (3.2)
Отметим, что эта величина не является
полноценной характеристикой гладкости, она
дает сигнал лишь о явном (существенном)
отклонении от гладкости. Эта величина равна
нулю, например, при изменении функции в
виде лесенки.
Для контроля по вычисленным значениям потенциала проверялась погрешность выполнимости теоремы Гаусса, согласно которой должно быть выполнено равенство
 En  ds  Q1  Q 2.
3. Обработка результатов и вычисление интегральных характеристик
(3.3)
S
Относительное отклонение от выполнимости этого равенства оценивалось по
формуле
dG  100(G1  G2  Q1  Q2) / 2/ Q1, (3.4)
Опишем вначале используемые подпрограммы. Итерации в стандартных узлах методом ПВР осуществлялись с использованием
п/п CalcV. При итерациях в нестандартных
узлах около сфер использовались п/п
в которой G1, G2 – вычисленные значения
потоков напряженности поля: G1 – поток че-
21
Е. Л. Тарунин
рез торцевые поверхности, G2 – поток через
боковую поверхность рассматриваемого объема. Эти потоки (интегралы) вычислялись по
формуле трапеций с аппроксимацией нормальных компонент напряженности центральными разностями на фиктивных границах, удаленных от внешних границ на 20 узлов пространственной сетки. В типичных расчетах
выполнялось
неравенство
G 2  3.5  G1 ; с ростом L и коэффициента
 это отношение увеличивалось. Величина
dG зависит от погрешностей вычисления потенциала, потоков G1, G2 и зарядов Q1,Q2.
Вычислительные эксперименты позволили выяснить, что основная погрешность dG
обусловлена погрешностью вычисления радиальных компонент напряженности поля на сферах En1[j], En2[j]. Частично этот вывод подтверждался негладкой зависимостью этих
функций от полярного угла. Для уменьшения
этого источника погрешности было испытано
несколько вариантов аппроксимации En1[j],
En2[j]. В итоге было выяснено, что наименьшая
погрешность вычисления En1[j], En2[j] соответствует варианту, в котором значения En1[j],
En2[j] вычислялись на поверхности сфер, удаленных от поверхности шаров на расстояние
ds=kh, которое немного больше шага пространственной сетки h (k>1). Значения потенциала на
этих сферах вычислялись по аппроксимации
значений потенциалов в ближайших точках
стандартной ячейки. В плоскости r, x в слой с
радиусами от R1 до (R1+kh) попадает около
2 k  RN стандартных ячеек с шагом h.
Первые подробные вычисления с различными значениями ds= k  h были выполнены при L  3,   10, RN  10. Минимальное значение характеристики отклонения
от гладкости GladEn =0.0678% было достигнуто при k=1.15 (при k=0 величина GladEn
=0.979%, а при k=1 GladEn=0.117%).
Был испытан также вариант вычисления
величин En1[j], En2[j] на трех радиусах
При вычислении потенциала на радиусе
(R1 +kh) была использована коррекция погрешности вычисления односторонней разностью величин En2[j], En2[j]. Эта коррекция
осуществлялась умножением полученных
значений En1[j], En2[j] на множитель
  1  ds / R . Множитель  вычислялся в
предположении, что потенциал вблизи сфер
близок к потенциалу для уединенной заряженной сферы. Описываемая коррекция снизила
( L  3,   10, RN  10 ) отклонение от выполнимости теоремы Гаусса до dG=0.21%, что
примерно в 55 раз меньше dG без коррекции.
Вариант расчетов потенциала на радиусе (R1+kh) позволил ликвидировать негладкость зависимостей En1[j], En2[j] и, самое
главное, снизить почти на порядок показатель
отклонения от выполнимости теоремы Гаусса
(типичные значения dG стали менее одного
процента для зарядов одного знака, симметрия решения задачи с зарядами разных знаков
при Q1=-Q2 приводила к dG=0).
Опишем формулы вычисления углов,
площадей, напряженности поля, зарядов и
горизонтальной компоненты силы для правой
половины первой сферы, полагая, что координаты j-го узла равны k, dm. Аналогичные
формулы используются для левых половин
сфер первого и второго шара. Индекс j введен
для нумерации нестандартных узлов, он изменяется от нуля до 2*JN1 (JN1=RN+p1). При
k<RN j=k, при k=RN учитываются нестандартные узлы для максимального значения k
на сфере. Номер j=JN1 соответствует углу
   / 2 , а номер j=2*JN1 – углу    .
Значения полярного угла вычислялись по
формуле
Mb1[ j ]  ArcTan(k / dm). (3.5)
Полагая, что границы соответствующего элемента площади заключены между углами b1,b2:
b1  0.5 Mb1 j   Mb1 j  1,
b 2  0.5 Mb1 j   Mb1 j  1,
ri  1  (k  0.025(i  2)) h (i=1,2,3)
3.6
вычисляются соответствующие элементы
площади (используется формула для шарового сегмента без площади основания)
с последующим усреднением с равными весами = 1/3. Характеристика гладкости при
этом практически не изменилась. Была надежда, что этот способ уменьшит амплитуду
колебаний kf при RN>24 (см. рис. 5). Эта надежда не оправдалась, так как полученные
значения kf изменились менее чем на
0.0017%.
Mds1[ j]  2* * R12 (cos(b1)  cos(b2)).(3.7)
Расстояние от j-го узла до поверхности
сферы
22
drj  h k 2  dm2  R1 позволяет
Задача электростатики о взаимодействие заряженных шаров …
L=3, вычисленное обычным способом, отличается от
приближенно вычислить радиальную компоненту поля
MEn1[ j ]  
Vs1  V [ k , mk 2[ j ]
. (3.8)
drj
FU  
на величину менее 0.01%.
В попытках уточнения решения испытывалось использование не нулевых значений
потенциала на внешней границе, равных
Характеристиками
зависимости
MEn1[j] являются среднее значение SEn и перепад напряженности kE1:
2*JN 1
SEn 
Q1 (3.2)  Q1 (2.8)
 Vs1
0.4
 MEn1[ j ]/(2 * JN1),
VG 
j 0
kE1  max MEn1[ j ]/ min MEn1[ j ]. (3.9)
1 Vs1 Vs 2
(

).
4 r1
r2
В этой формуле
По вычисленным значениям площадей
и напряженности поля вычисляются заряды
на элементах поверхности
Mdq1[ j ]  Mds1[ j ]* MEn1[ j ]. (3.10)
Суммирование значений Mdq1[j] позволяет вычислить заряды на обеих половинках
первого шара слева Q1L и справа Q1R.
Горизонтальная компонента силы, действующая на правую половину шара, вычисляется по формуле
r1
(3.12)
– расстояние от цен-
тра первого шара до границы, а r2 – расстояние от центра второго шара до границы. Использование значений (3.12) на границе
внешней области изменило коэффициент kf
примерно на 0.03 % при RN=10 и на 0.006%
при RN=16.
Малые изменения коэффициента kf были обнаружены и при использовании корректировки граничных значений с помощью
формулы экстраполяции
Vk ,0  2Vk ,1  Vk ,2 ,
JN1
f 12  0.5 Mdq1[ j]  MEn1[ j ]  cos(Mb1[ j ]).
которая следует из анализа поведения потенциала у границы. Последняя формула написана для левой границы области (x=0). Аналогичные формулы экстраполяции были использованы для других границ области.
j 0
(3.11)
Полная сила, действующая на первый
шар, равна F1=f11+f12. Эта сила сравнивается с силой, вычисляемой по закону Кулона
F 0  Q1 Q2/(4 L2 ) , с помощью коэффициента kf=F0/F. Значения kf>1 соответствуют
ослаблению силы по сравнению с законом
Кулона для точечных зарядов, а значения kf<1
– увеличению силы.
Часть формул, приведенных выше, относится к значениям номера j внутри соответствующего интервала; при значениях j=0,
j=JN1, j=2JN1 формулы корректируются.
Также меняются формулы и в случае, когда
значения напряженности поля вычисляются
на фиксированном радиусе r  R1  ds через
равные интервалы по углу    / 2 JN1,
Кроме вычисления сил по формулам
через значения нормальной компоненты поля
в ряде случаев использовалось вычисление
сил по значению производной от потенциальной энергии FU  U / L . Потенциальная
энергия одного шара вычислялась по формуле
U  Q1  Vs1 . Расчеты при L  3, 3  0.2 позволили выяснить, что значение силы при
Анализ графической зависимости радиальной компоненты поля En от угла  позволил обнаружить пульсации в зависимости
En(  ) . Была попытка сгладить эти пульсации, привлекая к вычислениям соседние значения по формуле
 ( j )  (1   )  En( j )    ( En( j  1)  En( j  1))
En
(3.13)
в которой  – весовой параметр, а индекс j
нумерует значения элементов массива
En( j ) . Наибольшее значение параметра
 =0.5 (при  =0 нет усреднения).
Сглаживание (3.1) эффективно при единичных пульсациях, а в целом нет, и поэтому
от него пришлось отказаться.
4. Результаты решения
Отметим, что более подробные данные
о результатах расчета содержатся в статье
[10]. Здесь же в основном приведены данные,
иллюстрирующие особенности решения.
23
Е. Л. Тарунин
На рис. 3 представлены изолинии потенциала для случая расстояния между центрами шаров L=3 при одинаковых радиусах
шаров и одинаковых потенциалах. Шаг между
изолиниями равнялся одной восьмой от максимального значения. Изображена не вся область решения (левый и правый край отстоят
от центров шаров на расстоянии, равном трем
радиусам). Поэтому в изображенную область
не вошли две изотермы.
По картине изолиний отчетливо видна
симметрия решения относительно середины
области. Взаимодействие шаров ослабило напряженность поля в пространстве между шарами. Этот эффект более отчетливо заметен
на зависимости радиальной компоненты поля
от угла (см. рис. 4). На этом рисунке изображена зависимость En(  ) для левого шара.
Массовым расчетам при произвольных
значениях расстояния между шарами предшествовали тестовые испытания при L=3. Тестовые испытания позволили определить параметры метода RN ,  , с помощью которых
можно получать достаточно точные результаты. Анализ этих данных дает возможность
утверждать, что погрешность величины kf
менее 1% достигается лишь при   20.
Серия расчетов при фиксированном
значении  =24 и различных значениях RN
позволила установить, что при RN  24:
1) модуль коэффициента kf колеблется
от 1.17697 до 1.17712 (иллюстрация этих колебаний представлена на рис. 5);
2) перепад напряженности kE1 колеблется около значения 2.8622 с относительной
амплитудой около 0.02%;
3) значение суммы потенциалов в граничных точках монотонно растет при увеличении RN, достигая при RN=40 значения
0.08075;
4) величина заряда монотонно растет
при увеличении RN, достигая при RN=40 значения 1.05327;
5) величины GladEn (характеристика
гладкости) и dG (погрешность отклонения от
выполнимости теоремы Гаусса) монотонно
убывают при увеличении RN.
Видно, что при малых углах (    / 2) напряженность значительно меньше, чем при
   / 2.
Для сравнения на рисунке изображена
горизонтальная линия со значением напряженности, равной напряженности для уединенного шара V=Vs1/(R1+h).
Заметим, что изображенная зависимость
в силу симметрии решения совпадает с зависимостью En(    ) для второго шара.
Большие значения

и RN приводят к
большим затратам машинного времени. Для
оценки машинного времени на РС с тактовой
частотой 2.1 GHz при   20 годится формула
t PC  1.035  10 5   2  RN 3 минут. (4.1)
Из этой формулы следует, что при
Рис. 3. Изолинии потенциала при L=3,
R2=R1, Vs1=Vs2
 =30 и
RN=30 время счета составляет 4.2 часа.
Зависимость коэффициента kf от числа
интервалов на радиусе RN представлена на
рис. 5. Как видно, при RN  20 значения kf
стабилизируются около 1.177 (это значение
показано на рисунке штриховой линией) с
отклонениями от него с относительной амплитудой менее 0.02%.
Рис. 4. Радиальная компонента напряженности поля
24
Задача электростатики о взаимодействие заряженных шаров …
1.178
Найденные уточненные значения kf0
(при уточнении счет выполнялся со значениями параметров метода RN ,  >24) представлены в таблице:
kf
1.177
1.176
1.175
L
2.25
2.5
2.75
kf0 1.44783 1.32195 1.23578
1.174
RN
1.173
0
10
20
30
L
3.0
3.5
4.0
kf0 1.17712 1.10687 1.06968
40
Рис. 5. Зависимость kf от числа интервалов
на радиусе RN (L=3,  =24)
Перейдем к обсуждению зависимости
коэффициента отношения сил kf+ от расстояния между центрами сфер для случая одноименных зарядов (Q1*Q2>0). Для нахождения
аналитической зависимости с двумя параметрами a0, c в виде
Анализ результатов, полученных при
различных значениях параметров метода RN
и  , позволил выяснить, что наибольшая погрешность решения (в основном по величине
коэффициента kf) обязана значениям параметра   30 . В справедливости этого утверждения можно убедиться путем анализа результатов для L=3, RN=20 и различных значений  .
Из анализа этих результатов следует зависимость
 kf  1.1756(1  0.0306 /  ).
kf   1  a0 * exp(c * ( L  2)),
(4.3)
использовались 10 значений (отобранные значения были получены для параметров метода
RN ,   30) :
L
kf
(4.2)
L
kf
Отсюда для получения значений kf с погрешностью менее 0.1% требуются значения
  31.
Из итоговых результатов тестирования
при L=3 следует, что достаточно точные результаты могут быть получены при значениях
 и RN  24. Колебания, представленные на
рис. 5, позволяют сделать вывод, что не всегда более высокие значения RN гарантируют
более точные результаты.
Часть вычислений kf при L, отличных от
3, проводилась следующим образом. При
фиксированном значении RN=24 выполнялся
счет для трех значений  (i ) =30, 25, 20 (i=1,
2, 3).
Затем по экстраполяции линейной зависимости kf  f 0  a  zi ( zi  1/  i ) на зна-
2.0
1.6273
2.75
1.2357
2.1
1.5478
3.0
117712
2.2
1.4789
2.25
1.4478
3.25
1.13594
3.5
1.10687
2.5
1.3219
4.0
1.06913
Функция (4.3) удовлетворяет асимптотическому стремлению kf  1 при L   .
Значения параметров этой функции a0, c находились методом градиентного спуска для
минимизации суммы квадратов невязок. Начальные значения параметров равнялись 1, а
шаги для их изменения da. dc равнялись
0.00001. Значения параметров аппроксимации
a0  0.64627, c  1.3239 , найденные после
353700 шагов градиентного спуска, снизили
начальную сумму квадратов невязок примерно в 120 раз до значения S=0.00553. Зависимость (4.3) изображена на рис. 6 штриховой
линией, крестиками отмечены табличные значения. Как видно, аппроксимация хороша для
значений L < 3.25, при больших значениях L
функция дает заниженные значения (наибольшее отклонение равно -2.2% при L=4).
По сумме квадратов невязок, и деленной на
число использованных табличных значений,
следует, что среднее отклонение табличных
данных от зависимости (4.3) около 0.4%.
Кроме функции (4.3) испытывалась аппроксимация также с двумя параметрами в виде
чение z=0 (   ) вычислялись два значения
kf 0 (1,2), kf 0 (1,3) , полученные для указанных
в скобках номеров расчетов i.
Далее полагалось, что
kf 0  (2 / 3)  kf 0 (1, 2)  (1/ 3)  kf 0 (1,3).
Оценка относительной погрешности вычислялась по формуле
kf  ( L)  1  a1/ L  a 2 / L2 .
  100(kf 0  kf1 ) / kf1.
25
(4.4)
Е. Л. Тарунин
ственной поправки. При уменьшении L эта
поправка возрастает.
3. Вычисленные поправки к закону
Кулона согласуются с расчетами других авторов, полученными другими методами, при
L=3; при L<3 поправочные коэффициенты
завышены, а при L>3 занижены.
4. Найдены приближенные аналитические зависимости для коэффициента поправки
к закону Кулона для зарядов одного знака
(для зарядов разного знака результаты представлены в работе [10]).
Эта функция имеет нужную асимптотику при L   . Однако получаемые значения коэффициентов имеют разные знаки, что
приводит на некотором интервале L к недопустимым значениям kf+  1 для зарядов с
одинаковыми знаками. Поэтому от этой аппроксимации пришлось отказаться.
1.8
kf
1.6
1.4
Список литературы
1.2
(L-2)
1. Саранин В.А. О взаимодействии двух
электрически заряженных шаров // Успехи
физических наук. 1999. Т. 169. С. 453–458.
2. Саранин В.А., Мейер В.В. Теоретические
и
экспериментальные исследования
взаимодействия двух проводящих заряженных шаров // Успехи физических наук,
2010. Т. 180, №10. С. 1009–1117.
3. Saranin V.A. Energy, force and field
strength in a system of two charged conducting balls // Journal of Electrostatics.
2013. Vol. 71. Р. 746–753.
4. Smythe W.R. Static and Dynamics Electrisity.
New York: McCraw – Hill, 1950.
5. Тарунин Е.Л. Электростатическое взаимодействие заряженных проводников на
близких расстояниях // Вестник Пермского университета. Сер. Физика. 2013.
Вып.2(24). С. 49–56.
6. Davis M.H. Two charged spherical conductors in a uniform electric field: forces and field
strength // Q.J. Mech. Appl. Math. 1964. Vol.
17. Р. 499–511.
7. Самарский А.А. Теория разностных схем.
М.: Наука, 1977. 656 с.
8. Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск, 1990. 226 с.
9. Физический энциклопедический словарь.
М.: Советская энциклопедия, 1960. Т. 5.
10. Тарунин Е.Л. Особенности электростатического взаимодействия заряженных сфер
на близких расстояниях // Вестник Пермского университета. Сер. Физика. 2014 (в
печати).
1
0
0.4
0.8
1.2
1.6
2
Рис. 6. Зависимость отношения сил от расстояния между сферами
Важной характеристикой зависимости
(4.3) является значение коэффициента kf+ при
dL  L  ( R1  R 2)  0 ,
соответствующее случаю соприкосновения
шаров. Из зависимости (4.3) следует, что это
предельное значение равно 1+a0=1.6462.
Счет по программе при L=2 дал чуть меньшее
значение kf+(2)=1.6245 0.0005. Линейная
экстраполяция на значение
dL  L  ( R1  R 2)  0 ,
построенная по двум значениям L=2.1 и
L=2.2, дает значение kf+=1.6157. Три различных способа вычисления предельного значения kf+(2) позволяют дать оценку этой величины kf+(2)  1.6216 с относительной погрешностью менее 0.4%.
Выводы
1. Реализован метод расчета электростатического поля около двух заряженных
шаров на основе решения уравнения Лапласа
для потенциала поля. Метод дает полную информацию о распределении потенциала и его
интегральных характеристиках.
2. Расчет показал, что на расстояниях
между центрами сфер L  2( R1  R2 ) закон
Кулона для точечных зарядов требует суще-
26
Задача электростатики о взаимодействие заряженных шаров …
Problem of electrostatic of the interaction
of two charged spheres at close distances
E. L. Tarunin
Perm State University, Russia, 614990, Perm, Bukireva st., 15
(342) 2 396-409
The electrostatic interaction between two conducting balls was studied by numerical method.
The strength of the interaction was found from the potential distribution of the electrostatic
field, which is determined by solving the finite-difference equations for the Laplace’s equation
in cylindrical coordinates. Used method to obtain complete information about the potential distribution described in detail. Analytical functions were built for describing the amendment to the
Coulomb’s law in the case of the same charges.
Key words: electrostatic; the Laplace’s equation; potential distribution of the electrostatic
field; the strength of the interaction conducting balls.
27
1/--страниц
Пожаловаться на содержимое документа