Технология расчета электромагнитных полей численными методами

Лабораторная работа № 5. Технология расчета электростатического поля силового кабеля методом конечных разностей

1. Цель работы

Ознакомление с методом конечных разностей на примере моделирования электростатического поля в силовом кабеле.

2. Теоретическая справка

Алгоритм задач электромагнитного поля методом конечных разностей:

1.На область задачи накладывается прямоугольная сетка. Область задачи состоит из точек и становится дискретной (рис. 1).

Рис. 1. Область задачи

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

3.Дифференциальные уравнения Пуассона и Лапласа с частными производными заменяются на разностные уравнения.

		−
2
	k		k −1

		k
k +1	−

		k +1			−
k	−		k	+		k

k ,
−1		k +1	− 2	k

+

	k −1

,

x x	k +1	− x	k

h

, y yk +1 − yk		h , где h – шаг дискретизации,
x2	h2	, y2 h2
h
	i,j+1
h
i-1,j	i,j	i+1,j

i,j-1

Рис. 2 Шаг дискретизации и координаты точек

Для каждой точки (i,j) внутри области формируется уравнение:

i+1, j

− 2

i, j

+

i−1, j

+

i, j+1

− 2

i, j

+

i, j−1

−

/ .

h

2

h

2

i, j

r

0

Для каждой граничной точки (i,j) с заданным потенциалом (условие Дирихле):

i, j гр .

Для каждой граничной точки (i,j) с n 0 (условие Неймана) на горизонтальной

оси:

i, j i+1, j ,

27

на вертикальной оси:

i, j

i,

j+1

.

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

Пример

Расчет электростатического поля заряженной двухпроводной линии. Плоскопараллельную задачу решаем для одной плоскости. Для любой другой параллельной плоскости решение будет совпадать.

y

U	=100В	U	=-100В
1		2

x

Рис. 3. Заряженная двухпроводная линия Нанесем на открытую область задачи сетку, искусственно ограничим ее границей. В

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

условию Неймана

n

0

.

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

	y				граничные точки с нулевым
	5					потенциалом
	4					внутренние точки
j	3
	2
	1					x
	1	2	3	4	5	6
				i
граничная точка с						граничные точки с
потенциалом U		=-100В				условием n 0
	2

Рис. 4. Расчетная область задачи

Запишем для каждой внутренней точки разностное уравнение (для снижения числа отображаемых переменных границу с нулевым потенциалом в уравнениях не учитываем).

(2,1)			2,1		2,2	,			Решение
(4,1)			4,1		4,2	,
(5,1)	5,1			5,2 ,
(2,2)	3,2			− 2 2,2			+	2,3 − 2 2,2 + 2,1	0 ,
								h2
				h2
									28

(3,2)

(4,2)

(5,2)

(2,3)

(3,3)

(4,3)

(5,3)

(2,4)

(3,4)

(4,4)

(5,4)

	4,2	− 2					3,2		+			2,2				+					3,3			− 2					3,2				− 100					0	,
						h	2																					h	2
	5,2	− 2					4,2		+		3,2					+					4,3			− 2					4,2				+			4,1		0	,
						h	2																					h	2
− 2 5,2						+ 4,2				+				5,3				− 2 5,2										+ 5,1						0 ,
			h2																			h2
	3,3	− 2					2,3		+			2,4				− 2						2,3				+			2,2				0 ,
		h		2																h		2
	4,3	− 2					3,3		+		2,3					+				3,4			− 2						3,3				+		3,2			0	,
						h	2																					h	2
5,3 − 2 4,3									+ 3,3							+		4,4					− 2 4,3										+ 4,2					0 ,
						h2																						h2
− 2			5,3			+		4,3		+					5,4			− 2						5,3				+			5,2			0			,
			h		2																	h		2
	3,4	− 2					2,4		+	− 2							2,4				+					2,3		0					,
		h		2													h		2
	4,4	− 2					3,4		+			2,4				+		− 2								3,4		+				3,3		0			,
						h	2																			h	2
	5,4	− 2					4,4		+		1,4					+		− 2							4,4			+				4,3		0			,
						h	2																			h	2
− 2 5,4						+ 4,4				+			− 2 5,4 + 5,3																	0 .
			h2
																			h2

Решение на сетке (серым отмечена расчетная область):

		2,1
		4,1
		5,1
		2,2
		3,2
		4,2
		5,2
		2,3
		3,3
		4,3
		5,3
		2,4
		3,4
		4,4
		5,4

			−14.5
			−15.5
			− 5.7
			−14.5
		− 35.9
			−15.5
			− 5.7
			− 7.8
			−13.4
			− 5.0
			−1.7
			− 3.1
			− 4.8
			− 2.7
			−1.1

В

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

—метода простых итераций — состоит в выполнении следующих процедур.

1.Исходная задача Ax = b преобразуется к равносильному виду:

x= αx+β

где α — квадратная матрица порядка n; β — столбец. Это преобразование может быть выполнено различными путями, но для обеспечения сходимости итераций нужно добиться выполнения условия ‖ ‖ .

2. Столбец β принимается в качестве начального приближения x(0) = β и далее многократно выполняются действия по уточнению решения, согласно рекуррентному соотношению

29

x(k+1) = αx(k)+β, k=0,1,2,…

или в развернутом виде

(	)	(	)	(	)	(	)
(	)	(	)	(	)	(	)
(	)	(	)	(	)	(	)

3. Итерации прерываются при выполнении условия

‖ ( )

( )‖

где ε >0 – заданная точность, которую необходимо достигнуть при решении задачи.

Пример

Получим систему линейных уравнений из предыдущего примера в виде x = αx+β:

(2,1)

(4,1)

(5,1)

(2,2)

(3,2)

(4,2)

(5,2)

(2,3)

(3,3)

(4,3)

(5,3)

(2,4)

(3,4)

(4,4)

(5,4)

(k +1) 2,1

(k +1) 4,1

(k +1) 5,1

(k +1) 2,2

(k +1) 3,2

(k +1) 4,2

5(,k2+1)

(k +1) 2,3

3(,k3+1)

(k +1) 4,3

(k +1) 5,3

(k +1) 2,4

3(,k4+1)

4(k,4+1)

5(,k4+1)

(k ) 2,2

(k ) 4,2

(k ) 5,2

1

4

1

4

1

4

14

1

4

14

1

4

1

4

1

4

141414

,

,

,

(k ) 3,2

(k ) 4,2

(k ) 5,2

(k ) 4,2

(k ) 3,3

(k ) 4,3

(k ) 5,3

(k ) 4,3

(k ) 3,4

(k ) 4,4

(k ) 5,4

(k ) 4,4

Начальное приближение

+	(k )	+		(k )		,
	3,2			2,1
+	(k )	+		(k )	+ 100				,
	2,2			3,3
+	(k )	+		(k )	+		(k )			,
	3,2			4,3			4,1
+	(k )	+		(k )	,
	5,3			5,1
+	(k )	+		(k )		,
	2,4			2,2
+	(k )	+		(k )	+		(k )	,
	2,3			3,4			3,2
+	(k )	+		(k )	+		(k )			,
	3,3			4,4			4,2
+	(k )	+		(k )		,
	5,4			5,2
+	(k )		,
	2,3
+	(k )	+		(k )	,
	2,4			3,3
+	(k )	+		(k )	,
	1,4			4,3
+	(k )	.
	5,3

3(,02) =100В, остальные потенциалы равны нулю.

30

Условия задачи

Рассмотрим силовой кабель на рабочее напряжение U с одной или двумя жилами (табл.1). В идеале проводники должны иметь общую ось или располагаться симметрично, однако в процессе изготовления или эксплуатации возможно смещение проводников на расстояние d, что может привести к пробою изоляции кабеля, если максимальная напряженность электростатического поля превысит допустимое значение.

3.Вопросы для коллоквиума

1.Каков алгоритм решения задач электромагнитного поля методом конечных разностей?

2.Как задаются граничные условия? Как задать однородное уравнение Неймана?

3.Опишите метод простой итерации? Как его применить для решения системы линейных алгебраических уравнений?

4.Назовите основные этапы расчета электростатического поля методом конечных разностей.

4.Рабочее задание

1.Нарисуйте в масштабе на бумаге в клетку (или на миллиметровке) модель согласно вариантам в табл. 1-4. Шаг дискретизации (шаг сетки) должен составлять не более 0,1R.

Таблица 1

	Номер бригады	Модель
	(номер компьютера)
		R	Жила c
		1	потенциалом U
			1

1, 4, 7, … 3n+1

R	Оболочка
	кабеля, U=0
	d

2, 5, 8, … 3n+2

	Жила 1 c
	потенциалом U
	1
R	Оболочка
	кабеля, U=0
R1	Жила 2 c
d	потенциалом -U	1

31

3, 6, 9, … 3n

			Жила 1 c
			потенциалом U
			1
R			Оболочка
			кабеля, U=0
d	R	1	Жила 2 c
			потенциалом -U
d			1

Таблица 2

	Номер группы		1, 2, 10, 18			3, 4, 11, 17			5, 6, 12, 16			7, 8, 15		9, 13, 14
	R, мм		10			20			30			40		45
	d		R/2			R/3			R/4			R/5		R/2
	Таблица 3
	Номер бригады	1, 6, 11, 16,			2, 7, 12, 17,			3, 8, 13, 18			4, 9, 14, 19		5, 10, 15, 20
	(номер компьютера)	21			22
	R1		d/2			d/2,5			d/3			d/3,5		d/4
	εr	2			2,2			2,4			2,6		3
	U1, В	550			500			750			380		220

2.Определите оси симметрии и антисимметрии, расчетную область задачи (см. лабораторную работу №1). Определить граничные и внутренние точки, граничные условия. Отсчет точек рекомендуется начать с левого верхнего угла. Тогда координаты точек будут совпадать со строками и столбцами в таблицах (матрицах) со значениями потенциалов и напряженностей электрического поля. Координаты по оси x – i будут соответствовать столбцам матриц, координаты по оси y – j будут соответствовать строкам матриц.

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

лабораторной работы в блокнот и сохраните в файл «название файла латиницей».m в свою папку, откройте его в программе MATLAB. Шаблон:

%Самостоятельно сформируйте матрицу type с элементами:

%0 - внутренняя точка, 1 - граничная точка с заданным потенциалом (условие

%Дирихле), 2 - граничная точка с dfi/dn = 0(условием Неймана)

%Например:

%[1 1 1 1 1 1 1 1 1 1 1;

%2 0 0 0 1 1 1 1 1 1 1;

%2 0 0 0 0 0 1 1 1 1 1;

%1 1 0 0 0 0 0 1 1 1 1;

%1 1 1 0 0 0 0 0 1 1 1;

%1 1 1 0 0 0 0 0 1 1 1;

%1 1 1 0 0 0 0 0 0 1 1;

%1 1 0 0 0 0 0 0 0 1 1;

%2 0 0 0 0 0 0 0 0 0 1;

32

%2 0 0 0 0 0 0 0 0 0 1;

%1 1 1 1 1 1 1 1 1 1 1;]

%и матрицу U заданного потенциала (для элементов с type = 1)

%Например:

% [0	0	0 0 0 0 0 0 0 0 0;
% 0	0	0 0 0 0 0 0 0 0 0;
% 0	0	0 0 0 0 0 0 0 0 0;

%550 550 0 0 0 0 0 0 0 0 0;

%550 550 550 0 0 0 0 0 0 0 0;

%550 550 550 0 0 0 0 0 0 0 0;

%550 550 550 0 0 0 0 0 0 0 0;

%550 550 0 0 0 0 0 0 0 0 0;

% 0		0	0 0 0 0 0 0 0 0 0;
% 0		0	0 0 0 0 0 0 0 0 0;
%	0	0	0	0 0 0 0 0 0 0 0;]

%--------------------------------------------------------------------------

% Введите сюда исходные данные (в системе СИ)

%--------------------------------------------------------------------------

h = _________;

% Матрицы type и U можно задать вручную с помощью команды new

%в окне Workspace и двойного клика на переменной. Если матрица создана %вручную, эти строки следует закомментировать (поставить знак %).

type = _______; U = __________;

%Начальное приближение

Phi_new = U;

%Метод простой итерации epsilon = 1;

I = size(type,1); J = size(type,2);

while epsilon >= 0.001

%Если цикл while не завершает работу, прервите его, нажав ctrl+C %или ctrl+break. Проверьте уравнения.

Phi_old = Phi_new;

for i = 1:I for j = 1:J

% Расчет для каждой точки в зависимости от ее типа

switch type(i,j) %при type=1 расчет проводить не надо, поскольку

%Phi(i,j) = U(i,j) case 0

%--------------------------------------------------------------------------

%введите сюда уравнение для внутренней точки

% Phi_new(i,j) = f(Phi_old); %--------------------------------------------------------------------------

Phi_new(i,j) = _________________; case 2

%--------------------------------------------------------------------------

33

%введите сюда уравнение для граничной точки с условием %Неймана

% Phi_new(i,j) = f(Phi_old); %--------------------------------------------------------------------------

Phi_new(i,j) = _____________;

end end

end %--------------------------------------------------------------------------

%введите сюда расчет epsilon

%--------------------------------------------------------------------------

epsilon = max(max(abs(____________)));

end

Phi = Phi_new;

% Расчет напряженности электрического поля

[Exneg, Eyneg] = gradient(Phi,h,h); Ex = -Exneg;

Ey = -Eyneg;

E = zeros(I,J); Ea = zeros(I,J); for i = 1:I

for j = 1:J %--------------------------------------------------------------------------

%Введите выражение для расчета модуля E и угла Ea вектора %напряженности электрического поля

%--------------------------------------------------------------------------

E(i,j) = _______________; Ea(i,j) = ______________;

end end

% Построение картины поля hx = 0:h:(I-1)*h;

hy = 0:h:(J-1)*h;

[x,y] = meshgrid(hx,hy);

contour(hx,hy,Phi), hold on, quiver(hx,hy,Ex,Ey), hold off set(gca,'YDir','reverse')

Выполните следующие задачи:

a.Создайте переменные type с типом точки: 0 – внутренняя точка, 1 – граничная точка с заданным потенциалом (условие Дирихле), 2 – граничная точка с условием n 0 (условие Неймана); и U для точек c заданным потенциалом. Переменные просматриваются и создаются в окне Workspace.

b.Введите исходные данные (в системе СИ);

new

c. Введите уравнение для нахождения i, j для внутренней точки согласно методу простой итерации;

34

d.Введите уравнение для нахождения условием Неймана;

new

i, j

для граничной точки с

e. Введите расчет точности полученного решения

new

;

f.Составляющие векторов напряженности электрического поля по оси x и y находятся с помощью функции gradient1. Введите формулы для расчета модуля и угла векторов напряженности электрического поля в точках расчетной области.

g.Запустите программу с помощью клавиши F5. Если возникают ошибки, внимательно проверьте выполнения всех пунктов.

4.В результате работы программы проводится расчет потенциалов и напряженности электрического поля в точках расчетной области, строится картина поля, например, как на рис. 5: Необходимо сравнить полученное решение с решением в программе ELCUT или FlexPDE.

y, м

0											550
0.001											500
0.002											450
0.003											400
0.004											350
0.005											300
0.006											250
0.007											200
0.008											150
0.009											100
0.01	0.001	0.002	0.003	0.004	0.005	0.006	0.007	0.008	0.009	0.01	50
0
					x,	м

Рис. 5. Картина поля

5.Зарисуйте картину поля во всей области задачи с учетом выбора осей симметрии (антисимметрии) с обозначением расположения жил и оболочки кабеля.

6.Определите точку с максимальной напряженностью электрического поля и величину максимальной напряженности электрического поля.

7.С помощью уравнения Пуассона Рассчитайте электрический заряд, наведенный на жилах на единицу длины кабеля. Принять, что объемный заряд в объеме h·h·l неизменен, а l = 1м. Сравните с результатом расчета электрического заряда в программе ELCUT или FlexPDE.

8.Рассчитайте электрическую емкость кабеля (между жилой и оболочкой или между жилами в зависимости от варианта) Сравните с результатом расчета электрической емкости в программе ELCUT или FlexPDE.

5. Вопросы к защите

1 Напряженность электрического поля связана с потенциалом выражением E = – gradφ.

35