ЗАДАЧА ЭЛЕКТРОСТАТИКИ О ВЗАИМОДЕЙСТВИИ ЗАРЯЖЕННЫХ ШАРОВ НА БЛИЗКИХ РАССТОЯНИЯХ Текст научной статьи по специальности «Физика»

2014

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

Математика. Механика. Информатика

Вып. 3(26)

УДК 537.2

Задача электростатики о взаимодействии заряженных шаров на близких расстояниях

Е. Л. Тарунин

Пермский государственный национальный исследовательский университет Россия, 614990, Пермь, ул. Букирева, 15 (342) 239-64-09

Известно, что сила взаимодействия заряженных проводников на близких расстояниях существенно отличается от закона Кулона для точечных зарядов, размещенных в центре шаров. Это отличие для проводящих сфер было подробно исследовано в работах [1-3], для тел цилиндрической формы - в [5]. В работах [1-3] расчеты были выполнены с использованием емкостных коэффициентов, вычисляемых с помощью бесконечных рядов [4]. В [5] вычисления были проведены на основе решения методом сеток задачи Дирихле для потенциала электростатического поля. В работе [6] для вычисления сил взаимодействия использовано решение уравнения Лапласа в бисферических координатах. В данном исследовании расчеты выполнены для проводящих сфер, как и в [1-3], но с использованием вычисляемого потенциала электростатического поля. Они показали как соответствие с результатами работ [1-3], так и отличия.

Ключевые слова: электростатика; уравнение Лапласа; потенциал электростатического поля; взаимодействие заряженных тел на близких расстояниях.

1. Постановка задачи

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

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

Во втором варианте основной счет выполнялся с использованием уравнения Лапласа, записанного в цилиндрической системе координат. Для облегчения обработки резуль-

© Тарунин Е. Л., 2014

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

Геометрия расчетной области изображена на рис. 1. Расчеты выполнялись с учетом осевой симметрии. Поэтому на рисунке изображена лишь половина области.

Рис. 1. Геометрия расчетной области

Параметрами задачи являются расстояние между центрами сфер L=x2 - х1, радиусы сфер R1, R2 и заряды на сферах Q1, Q2. Расчеты выполнялись в безразмерных переменных, в качестве единицы расстояния был выбран радиус первого шара R1 =1 (полагалось, что R2 — ^). Постановка задачи предполагает,

что взаимодействующие заряды изолированы. Положение внешней границы, на которой задавалось нулевое значение потенциала, определялось параметром метода /Л >10:

Х1 = АО= / +R2, 02С= / +R2,

АB=CD= ¡U+R2. (1.1)

Уравнение Лапласа для потенциала в цилиндрических координатах X, Г с учетом симметрии имело вид

I =0. (1.2)

г дг дг дх

В стандартных узлах уравнение аппроксимировалось на квадратной сетке Аг = Ах = h. Соответствующая система уравнений решалась итерационным методом. Особая аппроксимация уравнения (1.2) использовалась на оси и в узлах с нестандартными шаблонами вблизи шаров. Описание этих аппроксимаций дано при описании метода.

При заданных значениях потенциалов на сферах Vs1 и Vs2 из решения задачи Дирихле находилось распределение потенциала. При равных значениях потенциалов на сферах Vs2=Vs1 решалась задача о взаимодействии одноименных зарядов. Задание Vs2=-Vs1 приводило к задаче о взаимодействии разноименных зарядов.

По вычисленным значениям потенциалов находились различные характеристики решения: сила электростатического взаимодействия ¥, зависимость радиальной компоненты поля на сферах Ег от полярного угла Р, отношение максимальной радиальной компоненты поля на сфере к минимальной компоненте kE1 и другие.

Основной характеристикой решения являлось отношение вычисленной силы к силе взаимодействия по закону Кулона

kf = FJ F . Величина к/ показывает отклонение силы взаимодействия от закона Кулона для точечных зарядов (при к/ > 1 сила взаимодействия меньше F0, а при к/ < 1 больше).

Постановка задачи позволяла также определить электроемкость шаров. Согласно электростатике [9] в системе проводников заряды на проводниках связаны с потенциалами на них соотношением

(1.3)

Здесь С.. - коэффициенты емкости .-го тела, а С1 / - коэффициент электростатической индукции между телами с индексами у (этот коэффициент определяет величину заряда ц.,

когда потенциал /-го тела равен (., а все остальные проводники заземлены). Для коэффициентов в формуле (1.3) выполняются соотношения

См > 0, С„/ = С/,. < 0 при . Ф / (1.4)

В случае двух проводников уравнения (1.3) имеют вид

41 = С1,1 ( + С1,2 (

42 = С2,1 ■ (1 + С2,2 ■ (2.

(1.5)

Коэффициенты в этих формулах зависят от радиусов сфер и от расстояния между сферами. При равенстве радиусов (именно этот

случай рассматривается в статье) С11 = С2 2 .

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

/ С1.1 + С1,2((2/(1) п „ а = 4Х142 =-1——-—. (16)

С1,2 + С2,2((2/(1)

Решение сформулированной задачи с потенциалами р1 = 1, р2 = 0 и вычисление

зарядов позволяют из формул (1.5) определить коэффициенты

С1,1 = Ц1, С1,2 = Ц2.

(1.7)

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

42, С1 2 стремятся (как и положено) к нулю, а

значения q^í, C^í у в выбранных единицах измерения - к 4к (при L=6, например, отличие q^í, С 1 от 4п составляет 8%). Вычисленные

коэффициенты позволяют найти отношение зарядов

У =

q + / q ~

C1,1 + С1.2 С1,1 — С1,2

< l. (18)

Здесь q соответствует случаю одноименных зарядов (( • (р2 = 1), а q - разноименных (( • (р2 =—1). При удалении зарядов у — 1, а c12 —0.

Основными параметрами сетки служили число интервалов на первом радиусе RN и параметр /Л, определяющий внешний размер

области (см. (1.1)). Максимальное значение индекса k равнялось NMR, а максимальное значение индекса m равнялось NM. Для хранения всех элементов массива для потенциала V[k,m]требовалось

n=RN2(L+2(R2+ / ))(R2+ Л )

ячеек оперативной памяти. При типичных значениях параметров задачи и метода (L=3, RN=30, / =24, R2=1) число элементов массива для потенциала электростатического поля n более миллиона.

2. Метод решения

-((* +л/2) • V -о; -л/2) • V) V =о. (2.1) Гк

Замена р ^ V сделана из методических соображений для отличия разностного решения от точного решения уравнения (1.2). После введения обозначений

БЩ =0.25(1+0.5/к), БЩ =0.25(1-0.5/к) (2.2)

уравнение (2.1) записывалось в виде, удобном для применения итерационного метода

Щт] =(В[к] • У[к+1,m] +B2[k] • У[к—1,m]+ У[к,т+1]+У[к,ш—1])/4.

(2.3)

Для сокращения объема вычислений значения компонент массивов В1 и В2 вычислялись один раз перед началом итераций.

Система уравнений для потенциала электростатического поля решалась методом последовательной верхней релаксации [7, 8]. Для узлов со стандартным шаблоном параметр верхней релаксации вычислялся по формуле, подобранной экспериментальным путем,

со,= 2/(1 + sm(0.9W NMR)). (2.4)

В этой формуле NMR - максимальное значение индекса k по радиусу.

Для узлов с нестандартным шаблоном вблизи сфер параметр релаксации равнялся 1 (или был чуть больше 1).

Первое слагаемое в уравнении Лапласа (1.2) на оси имеет устранимую неопределенность типа ноль деленный на ноль. Применение к этому слагаемому правила Лопиталя позволяет выяснить, что

lim

„0(14 (r ^)) = 2( К

r or or or

(2.5)

и использовать с учетом симметрии аппроксимацию

02( 4(VLm — У„ш )

2( ) -

(2.6)

Аппроксимация уравнения Лапласа (1.2) на квадратной сетке

rk = к • h, Хт = Ш • h в обозначениях школы

А.А. Самарского [7] имела вид

Полная аппроксимация уравнения Лапласа на оси с учетом (2.6) приобретает вид

^ = (4^ + ^ + ^-^/6. (2.7)

Для сокращения общего числа итераций

значения потенциала на оси

V

0,ш

находились

методом скалярной прогонки в предположении, что значения V т известны. Коэффициенты скалярной прогонки Рт, Qm вычислялись по рекуррентным формулам

1

6 — Pm—1

, Q = P(4V + Q 1) (2.8)

' ^m mV 1,m z^m—1/

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

W

N

Р

W е

S

к

Рис. 2. Обозначения для узла с нестандартным шаблоном

Выпишем конечно-разностную аппроксимацию уравнения Лапласа для нестандартного шаблона, изображенного на рис. 2.

Для обозначения узлов использована географическая символика с заглавными буквами ^Е^^ ^ - северный узел, W - западный, Е - восточный, S - южный; соответствующие маленькие буквы обозначают длину плеч в единицах стандартного шага

Для узлов, расположенных слева от сфер, 1, а для узлов расположенных справа, е=1. В использованных обозначениях аппроксимация уравнения Лапласа имела вид

-Г

rk

h

V-V

^y(Q%hn-s-h}+

s-h

(2.9)

1

(У-У У-У^

03e■h+w■h) &И м>Н

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

= k -у/ RN2 - dm2, = dm -^RN2 - k2.

(2.10)

Формулы приведены для узлов, расположенных справа от первой сферы. Подобный вид имеют формулы для второй сферы и для узлов, расположенных слева от сфер.

Наличие нестандартных узлов отражается и на организации циклов по узлам со стандартным шаблоном. В программе вычислялись компоненты четырех векторов тк1[к], тк2[к], тк3[к], тк4[к], которые определяли значение индекса по горизонтали т с нестан-

дартным шаблоном, наиболее удаленным от поверхности соответствующей сферы. Первые два вектора относятся к первой сфере, а следующие два - ко второй.

Итерации в стандартных узлах выполнялись с использованием подпрограммы с именем "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. Для нестандартного узла вблизи пер-

S

вого шара радиальная компонента поля вычисляется по формулам

<Г=-В1, Ег = (У[к,т]-Ш)/<Г. (2.11)

Опишем используемую классификацию нестандартных узлов. Назовем узлы с индексами тк1[к], тк2[к], тк3[к], тк4[к] для заданных значений к "граничными" нестандартными узлами (присвоение этим узлам названия граничных связано с тем, что рядом с ними слева или справа находится узел со стандартным шаблоном).

Анализ показывает, что для некоторых значений к близких к k=RN есть нестандартные узлы, которые расположены ближе к сфере. Такие узлы будем называть "близнецами". Число близнецов при фиксированном k<RN может быть оценено по формуле, в которой использована процедура вычисления целой части от вещественного числа:

тах <т = 1гипс(2^ М -1 -V 2 ■ М-1) .(2.12)

При Ш=20 тах йт =2, а при Ш=100 тах йт =5. Напомним, что число близнецов указывается для узлов, расположенных с одной стороны сферы (общее число близнецов вдвое больше).

Наибольшее число близнецов при k=RN. Близнецы при k=RN выделены в отдельную группу. Максимальное число близнецов при k=RN равно

р1 = 1птс(4 2 ■ т -1 + 0.000001). (2.13)

Отсюда следует, что при RN= 10 р1=4, а при RN=20 р1=6.

Перейдем к описанию основных блоков программы.

1. В первом блоке описываются используемые переменные, подпрограммы и задаются значения параметров задачи и метода.

2. Во втором блоке вычисляются параметры нестандартных узлов: индексы граничных нестандартных узлов йтк1[к], <тк2[к], показывающих удаление узла от оси сфер; длины плеч нестандартных узлов ^к1[к], hsk1[k] для первого шара и ^к2[к], hsk2[k] для второго; р1, р2 - число нестандартных узлов при k=RN и k=RN2 соответственно; Ы1[к[, к<2[к] - число близнецов. Заметим, что часть массивов, упомянутых в этом блоке, относятся лишь к правой стороне сфер. Для левых половин соответствующие значения вычисляются из соображений симметрии.

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

s1 := (B1[k] • V[k +1, m] + B2[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.

Заметим, что величина невязки nev вычислялась также в узлах с нестандартными шаблонами. При завершении 4-го блока счетчик внешних итераций mfull увеличивался на 1 и выполнялось сравнение невязки с eps0=10-13. Возврат в начало 4-го блока происходил при одном из двух условий -m0 < NM, nev>eps0. Первое условие экономило время счета, сокращая число итераций на первых внешних циклах. Без первого условия общее число внутренних итераций Smo возрастало примерно на 40%. Обычно до mfull ~ 40 число итераций m0=NM, а затем убывало до 1. Так, например, при

L = 3, ¡л = 10, RN = 45, eps = 1013

зависимость числа внутренних итераций то от числа внешних итераций т/иП указана в таблице:

mfull < 30 40 50 60 70 80 > 90

mo 1175 930 279 57 3 2 1

Обычно значения mo выводились с интервалом для внешних итераций, равном 10. Для интегрального анализа сходимости на печать выводились номера полных итераций ml, m2 (ш1 соответствует номеру, при котором впервые mo < NM, m2 соответствует номеру, после которого все последующие mo=1). Так, например, при значениях L=3, / =30, RN=24 полное число итераций

Smo=73 424, а номера ml, m2 равны соответственно 40 и 102. Полное число внутренних итераций Smo возрастало при увеличении числа интервалов на первом радиусе RN

сильнее линейной зависимости (Smo ~ RN12).

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

2JN1

SSE(mfull) = £ MEn1[ j ]. (2.14)

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] . По значениям радиальных компонент поля находятся отношения максимальной компоненты поля к минимальной kEl, kE2. Затем вычисляются заряды на половинках сфер QR1, QR2, QL1, QL2 и силы взаимодействия fll, f12, f21, f22. В конце обработки результатов расчета на экране строится зависимость радиальной компоненты поля от полярного угла. Завершается программа построением изолиний.

3. Обработка результатов и вычисление интегральных характеристик

Опишем вначале используемые подпрограммы. Итерации в стандартных узлах методом ПВР осуществлялись с использованием п/п CalcV. При итерациях в нестандартных узлах около сфер использовались п/п

CaseLeft, CaseRight; параметр верхней релаксации для узлов около сфер был близок к 1.

По вычисленным значениям потенциала V[k,m] находились следующие характеристики решения: радиальные компоненты поля на сферах Еп1[}], Еп2[}], заряды на полусферах QL1, QR1, QL2, QR2, полные заряды на сферах Q1=QL1+QR1, Q2=QL2+QR2, силы /11, /12/21/22, F1=/11+/12, F2=/21+/22. Значения Еп1[}], Еп2[}] позволяли вычислить как заряды на сферах, так и компоненты сил. Интегральными характеристиками зависимости Еп1[}] служили три величины - кЕ1 (перепад напряженности на сфере), SEn: (среднее значение) и GladEn (характеристика гладкости):

кЕ1( }) = тах(Еп1[ }])/тт(Еп1[ }]),

] }

SEn = —1— V Еп1[} ]. (3.1) 2 • JN1¿-J

Характеристика гладкости была введена для того, чтобы сравнивать алгоритмы сглаживания. В первых расчетах эффекты негладкой зависимости напряженности от полярного угла определялись по виду графического изображения функции Еп1(Р). Позднее была введена числовая характеристика GladEn. Для ее вычисления в цикле по индексу } от 2 до (2Ж1-1) суммировались абсолютные значения выхода Еп1[}] за пределы интервала от Еп1[}-1] до Еп1[}+1]. Затем в процентах вычислялось окончательное значение

GladEn := 100GladEn / 8Еп /(Л -1). (3.2)

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

Для контроля по вычисленным значениям потенциала проверялась погрешность выполнимости теоремы Гаусса, согласно которой должно быть выполнено равенство

[|] Еп • ds = 01 + Q 2. (3.3)

Относительное отклонение от выполнимости этого равенства оценивалось по формуле

dG = 100(01 + 02 - 01 - 02)/2/01, (3.4)

в которой G1, G2 - вычисленные значения потоков напряженности поля: G1 - поток че-

рез торцевые поверхности, G2 - поток через боковую поверхность рассматриваемого объема. Эти потоки (интегралы) вычислялись по формуле трапеций с аппроксимацией нормальных компонент напряженности центральными разностями на фиктивных границах, удаленных от внешних границ на 20 узлов пространственной сетки. В типичных расчетах выполнялось неравенство О2 > 3.5 • 01; с ростом L и коэффициента /Л это отношение увеличивалось. Величина dG зависит от погрешностей вычисления потенциала, потоков 01, 02 и зарядов 01,02.

Вычислительные эксперименты позволили выяснить, что основная погрешность dG обусловлена погрешностью вычисления радиальных компонент напряженности поля на сферах Еп1[}], Еп2[}]. Частично этот вывод подтверждался негладкой зависимостью этих функций от полярного угла. Для уменьшения этого источника погрешности было испытано несколько вариантов аппроксимации Еп1}], Еп2[}]. В итоге было выяснено, что наименьшая погрешность вычисления Еп1[}], Еп2[}] соответствует варианту, в котором значения Еп1[}], Еп2[}] вычислялись на поверхности сфер, удаленных от поверхности шаров на расстояние ds=kh, которое немного больше шага пространственной сетки Л (к>1). Значения потенциала на этих сферах вычислялись по аппроксимации значений потенциалов в ближайших точках стандартной ячейки. В плоскости г, х в слой с радиусами от R1 до (R1+кЛ) попадает около 2пк • RN стандартных ячеек с шагом Л.

Первые подробные вычисления с различными значениями ds= к • Л были выполнены при Ь = 3, /= 10, Ш = 10. Минимальное значение характеристики отклонения от гладкости 0^Еп =0.0678% было достигнуто при к=1.15 (при к=0 величина 0^Еп =0.979%, а при к=1 0ШЕп=0.117%).

Был испытан также вариант вычисления величин Еп1[}], Еп2[}] на трех радиусах

Г = 1 + (к + 0.025(/ - 2))Л (1=1,2,3)

с последующим усреднением с равными весами = 1/3. Характеристика гладкости при этом практически не изменилась. Была надежда, что этот способ уменьшит амплитуду колебаний к/ при RN>24 (см. рис. 5). Эта надежда не оправдалась, так как полученные значения к/ изменились менее чем на 0.0017%.

При вычислении потенциала на радиусе (R1 +кЛ) была использована коррекция погрешности вычисления односторонней разностью величин Еп2[}], Еп2[}]. Эта коррекция осуществлялась умножением полученных значений Еп1[}], Еп2[}] на множитель р = 1 + ds / R . Множитель р вычислялся в

предположении, что потенциал вблизи сфер близок к потенциалу для уединенной заряженной сферы. Описываемая коррекция снизила (Ь = 3, / = 10, RN = 10) отклонение от выполнимости теоремы Гаусса до dG=0.21%, что примерно в 55 раз меньше dG без коррекции.

Вариант расчетов потенциала на радиусе (Rl+kh) позволил ликвидировать негладкость зависимостей Еп1[}], Еп2[}] и, самое главное, снизить почти на порядок показатель отклонения от выполнимости теоремы Гаусса (типичные значения dG стали менее одного процента для зарядов одного знака, симметрия решения задачи с зарядами разных знаков при Q1=-Q2 приводила к dG=0).

Опишем формулы вычисления углов, площадей, напряженности поля, зарядов и горизонтальной компоненты силы для правой половины первой сферы, полагая, что координаты }-го узла равны к, dm. Аналогичные формулы используются для левых половин сфер первого и второго шара. Индекс} введен для нумерации нестандартных узлов, он изменяется от нуля до 2*^1 (JN1=RN+p1). При k<RN }=к, при k=RN учитываются нестандартные узлы для максимального значения к на сфере. Номер j=JN1 соответствует углу (3 = я72, а номер }=2*Ж1 - углу Р = п . Значения полярного угла вычислялись по формуле

МЬ1[}] = АгсТап(к / dm). (3.5)

Полагая, что границы соответствующего элемента площади заключены между углами Ь1,Ь2:

Ъ1 = 0.5* (М&1[/] + МЪ1\}} -1]),

Ь2 = 0.5*(МЪ1\} ] + МЪ1[} +1]), (3.6)

вычисляются соответствующие элементы площади (используется формула для шарового сегмента без площади основания)

Mds1[j] = 2* п* R12(cos(b1) - со8(Ъ2)). (3.7)

Расстояние от }-го узла до поверхности

сферы = лук + dm - R1 позволяет

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

MEn1[ j ] = - Vs1 - V 'k •mk 2'j ]. (3.8)

drj

Характеристиками зависимости

MEn1[j] являются среднее значение SEn и перепад напряженности kE1:

2JN1

SEn = £ MEn1[ j] /(2 * JN1),

j=0

kE1 = max MEn1[ j ]/min MEn1[ j ].(3.9)

По вычисленным значениям площадей и напряженности поля вычисляются заряды на элементах поверхности

Mdq 1[ j ] = Mds 1[ j ] * MEn1[ j ]. (3.10)

Суммирование значений Mdq1£j] позволяет вычислить заряды на обеих половинках первого шара слева Q1L и справа Q1R.

Горизонтальная компонента силы, действующая на правую половину шара, вычисляется по формуле

JN1

f 12 = 0.5^Mdq1[j] • MEnl[ j] • cos(M>1[j ]).

j=0

(3.11)

Полная сила, действующая на первый шар, равна F1=f11+f12. Эта сила сравнивается с силой, вычисляемой по закону Кулона F 0 = Q1■ Q2/(4лL2), с помощью коэффициента kf=F0/F. Значения к/>1 соответствуют ослаблению силы по сравнению с законом Кулона для точечных зарядов, а значения к<1 - увеличению силы.

Часть формул, приведенных выше, относится к значениям номера у внутри соответствующего интервала; при значениях у=0, ]'=^1, ]=2М1 формулы корректируются. Также меняются формулы и в случае, когда значения напряженности поля вычисляются на фиксированном радиусе г = R1 + ds через равные интервалы по углу А/ = п /(2JNl),

Кроме вычисления сил по формулам через значения нормальной компоненты поля в ряде случаев использовалось вычисление сил по значению производной от потенциальной энергии ^ = -ди / дL . Потенциальная энергия одного шара вычислялась по формуле и = Q1 ■Vs1. Расчеты при L = 3, 3 ± 0.2 позволили выяснить, что значение силы при

L=3, вычисленное обычным способом, отличается от

=_ а(з.2) - а(2.8) _

и 0.4

на величину менее 0.01%.

В попытках уточнения решения испы-тывалось использование не нулевых значений потенциала на внешней границе, равных

1 Ш Vs2

Уо =— (— +-).

4л г Г

(3.12)

В этой формуле Г1 - расстояние от центра первого шара до границы, а Г - расстояние от центра второго шара до границы. Использование значений (3.12) на границе внешней области изменило коэффициент к/ примерно на 0.03 % при Ш=10 и на 0.006% при RN=16.

Малые изменения коэффициента к/ были обнаружены и при использовании корректировки граничных значений с помощью формулы экстраполяции

V ,0 = 2Ук ,1 - Ук ,2,

которая следует из анализа поведения потенциала у границы. Последняя формула написана для левой границы области (х=0). Аналогичные формулы экстраполяции были использованы для других границ области.

Анализ графической зависимости радиальной компоненты поля Еп от угла ¡3 позволил обнаружить пульсации в зависимости Еп(3). Была попытка сгладить эти пульсации, привлекая к вычислениям соседние значения по формуле

Еп(у) = (1 - 3) ■ Еп(у) + 3 ■ (Еп(у -1) + Еп(у +1))

(3.13)

в которой 3 - весовой параметр, а индекс у нумерует значения элементов массива Еп(у). Наибольшее значение параметра

3=0.5 (при 3=0 нет усреднения).

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

4. Результаты решения

Отметим, что более подробные данные о результатах расчета содержатся в статье [10]. Здесь же в основном приведены данные, иллюстрирующие особенности решения.

На рис. 3 представлены изолинии потенциала для случая расстояния между центрами шаров Ь=3 при одинаковых радиусах шаров и одинаковых потенциалах. Шаг между изолиниями равнялся одной восьмой от максимального значения. Изображена не вся область решения (левый и правый край отстоят от центров шаров на расстоянии, равном трем радиусам). Поэтому в изображенную область не вошли две изотермы.

По картине изолиний отчетливо видна симметрия решения относительно середины области. Взаимодействие шаров ослабило напряженность поля в пространстве между шарами. Этот эффект более отчетливо заметен на зависимости радиальной компоненты поля от угла (см. рис. 4). На этом рисунке изображена зависимость Еп(() для левого шара.

Видно, что при малых углах (( < п/2)на-пряженность значительно меньше, чем при Р ^ п/2.

Для сравнения на рисунке изображена горизонтальная линия со значением напряженности, равной напряженности для уединенного шара V=Vs1/(R1+h).

Заметим, что изображенная зависимость в силу симметрии решения совпадает с зависимостью Еп(п - () для второго шара.

Рис. 3. Изолинии потенциала при L=3, R2=R1, Vs1=Vs2

Рис. 4. Радиальная компонента напряженности поля

Массовым расчетам при произвольных значениях расстояния между шарами предшествовали тестовые испытания при L=3. Тестовые испытания позволили определить параметры метода RN, ¡1, с помощью которых

можно получать достаточно точные результаты. Анализ этих данных дает возможность утверждать, что погрешность величины 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.

Большие значения ¡1 и RN приводят к

большим затратам машинного времени. Для оценки машинного времени на РС с тактовой частотой 2.1 GHz при ¡1 > 20 годится формула

tPC «1.035 • 10 5 • ц2 • RN3 минут. (4.1)

Из этой формулы следует, что при ¡=30 и RN=30 время счета составляет 4.2 часа.

Зависимость коэффициента kf от числа интервалов на радиусе RN представлена на рис. 5. Как видно, при RN > 20 значения kf стабилизируются около 1.177 (это значение показано на рисунке штриховой линией) с отклонениями от него с относительной амплитудой менее 0.02%.

1.178 —

- И 1.177------

1.176 —

1.175 —

1.174 —

- RN 1.173--1-1-1-1-1-1-1-1

0 10 20 30 40

Рис. 5. Зависимость kf от числа интервалов на радиусе RN (L=3, ц =24)

Анализ результатов, полученных при различных значениях параметров метода КЫ и ц, позволил выяснить, что наибольшая погрешность решения (в основном по величине коэффициента к/) обязана значениям параметра ц < 30. В справедливости этого утверждения можно убедиться путем анализа результатов для Ь=3, КЫ=20 и различных значений ц.

Из анализа этих результатов следует зависимость

-к/ « 1.1756(1 + 0.0306/ц). (4.2)

Отсюда для получения значений к/ с погрешностью менее 0.1% требуются значения ц > 31.

Из итоговых результатов тестирования при Ь=3 следует, что достаточно точные результаты могут быть получены при значениях ц и КЫ> 24. Колебания, представленные на

рис. 5, позволяют сделать вывод, что не всегда более высокие значения КЫ гарантируют более точные результаты.

Часть вычислений к/ при Ь, отличных от 3, проводилась следующим образом. При фиксированном значении КЫ=24 выполнялся счет для трех значений ц(/ )=30, 25, 20 (/=1,

2, 3).

Затем по экстраполяции линейной зависимости к/ = /0 + а • zi (zi = 1/ ц ) на значение z=0 (ц = да) вычислялись два значения к/0 (1,2), к/0 (1,3) , полученные для указанных в скобках номеров расчетов /.

Далее полагалось, что

к/0 = (2/3) • /,(1,2) + (1/3) • £/,(1,3).

Оценка относительной погрешности вычислялась по формуле

3 = 100(/ - к/,)/к/г.

Найденные уточненные значения к/0 (при уточнении счет выполнялся со значениями параметров метода КЫ, ц >24) представлены в таблице:

L 2.25 2.5 2.75

к!о 1.44783 1.32195 1.23578

L 3.0 3.5 4.0

к!о 1.17712 1.10687 1.06968

Перейдем к обсуждению зависимости коэффициента отношения сил к/+ от расстояния между центрами сфер для случая одноименных зарядов (Ql*Q2>0). Для нахождения аналитической зависимости с двумя параметрами а0, с в виде

/ = 1 + а0* ехр(-с *(Ь - 2)), (4.3)

использовались 10 значений (отобранные значения были получены для параметров метода

т, ц> 30):

L 2.0 2.1 2.2 2.25 2.5

к! 1.6273 1.5478 1.4789 1.4478 1.3219

L 2.75 3.0 3.25 3.5 4.0

к:! 1.2357 117712 1.13594 1.10687 1.06913

Функция (4.3) удовлетворяет асимптотическому стремлению к/ ^ 1 при Ь ^ да .

Значения параметров этой функции а0, с находились методом градиентного спуска для минимизации суммы квадратов невязок. Начальные значения параметров равнялись 1, а шаги для их изменения da. dc равнялись 0.00001. Значения параметров аппроксимации а0« 0.64627, с «1.3239, найденные после 353700 шагов градиентного спуска, снизили начальную сумму квадратов невязок примерно в 120 раз до значения S=0.00553. Зависимость (4.3) изображена на рис. 6 штриховой линией, крестиками отмечены табличные значения. Как видно, аппроксимация хороша для значений Ь < 3.25, при бо'льших значениях Ь функция дает заниженные значения (наибольшее отклонение равно -2.2% при Ь=4). По сумме квадратов невязок, и деленной на число использованных табличных значений, следует, что среднее отклонение табличных данных от зависимости (4.3) около 0.4%. Кроме функции (4.3) испытывалась аппроксимация также с двумя параметрами в виде

/ (Ь) = 1 + а1/ Ь + а 2/Ь2. (4.4)

Эта функция имеет нужную асимптотику при Ь —> го . Однако получаемые значения коэффициентов имеют разные знаки, что приводит на некотором интервале Ь к недопустимым значениям к/+ < 1 для зарядов с одинаковыми знаками. Поэтому от этой аппроксимации пришлось отказаться.

Рис. 6. Зависимость отношения сил от расстояния между сферами

Важной характеристикой зависимости (4.3) является значение коэффициента к/+ при

dL = Ь - (R1 + R2) = 0,

соответствующее случаю соприкосновения шаров. Из зависимости (4.3) следует, что это предельное значение равно 1+а0=1.6462. Счет по программе при Ь=2 дал чуть меньшее значение к/+(2) = 1.6245 ±0.0005. Линейная экстраполяция на значение

dL = Ь - (R1 + R2) = 0,

построенная по двум значениям Ь=2.1 и Ь=2.2, дает значение к/+ = 1.6157. Три различных способа вычисления предельного значения к/+(2) позволяют дать оценку этой величины к/+(2) ~ 1.6216 с относительной погрешностью менее 0.4%.

Выводы

1. Реализован метод расчета электростатического поля около двух заряженных шаров на основе решения уравнения Лапласа для потенциала поля. Метод дает полную информацию о распределении потенциала и его интегральных характеристиках.

2. Расчет показал, что на расстояниях между центрами сфер Ь < 2(R + R2) закон Кулона для точечных зарядов требует суще-

ственной поправки. При уменьшении L эта поправка возрастает.

3. Вычисленные поправки к закону Кулона согласуются с расчетами других авторов, полученными другими методами, при L=3; при L<3 поправочные коэффициенты завышены, а при L>3 занижены.

4. Найдены приближенные аналитические зависимости для коэффициента поправки к закону Кулона для зарядов одного знака (для зарядов разного знака результаты представлены в работе [10]).

Список литературы

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 (в печати).

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.