1.Как согласно методу конечных разностей формируется уравнение Пуассона, Лапласа для определенной точки?

2.Чем внутренние точки отличаются от граничных?

3.Как рассчитать напряженность электрического поля в точке, если известен потенциал во всех точках рассматриваемой области?

4.Каково начальное приближение для расчета потенциалов было использовано в работе?

5.Как рассчитать электрическую емкость кабеля?

6.Как рассчитать электрический заряд?

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

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

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

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

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

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

21		22		23		24	25
e	25	e	26	e	27		e
							28
	e		e		e		e
	29		30		31		32
16		17		18		19	20
e	17	e	18	e	19		e
							20
	e		e		e		e
	21		22		23		24
		100		В
11		12		13		14	15

e		e	10		e		e		12
9					11
e13			e14		e15				e16
6		7		8		9			10
e		e	2		e		e	4
1					3
e	5		e		e	7			e
			6						8
1		2		3		4			5

Рис. 1. Область задачи, разбитая на конечные элементы

2.Формируется список конечных элементов в соответствии с обозначением вершин i, j, k на рис. 2. в виде

Конечный элемент	i	j	k

36

	y
y	, y	k	i	k
i

	h
y
j	j

			h
x	, x	j	x	k
i

y

yj

y , y	k
i

x

j

h

k				i
	h
xk	x	, x	j	x
	i

Рис. 2 Координаты а) верхнего и б) нижнего треугольников.

3.Для целей формирования матричного уравнения вида Ax=b, поэлементно для каждой строки из списка конечных элементов в матрицу коэффициентов Ap добавляются коэффициенты каждого из треугольников. Сначала формируется нулевая mxm матрица, где m – число точек (вершин треугольников):

A	(0)		0
	p

.

Далее на каждом шаге l добавляется вклад каждого из треугольников

i

j

k

i

1

− 0,5

− 0,5

A

(l )

A

(l −1)

− 0,5

0,5

0

+ j

.

p

p

− 0,5

0

0,5

k

или

A(l ) A(l −1) + 1;

ii

ii

A(l ) A(l −1)

− 0,5; A(l ) A(l −1)

− 0,5;

ij

ij

pi

ji

A(l ) A(l −1)

− 0,5; A(l ) A(l −1)

− 0,5;

ik

ik

ki

ki

A(l )

A(l −1)

+ 0,5;

jj

jj

A(l )

A(l −1)

+ 0,5.

kk

kk

37

4.Строки, соответствующие граничным точкам, в которых известен потенциал, вычеркиваются. Столбцы, соответствующие граничным точкам, домножаются на значения известных потенциалов и со знаком минус переносятся в правую часть уравнения. Для двух граничных точек s и r это будет выглядеть следующим образом:

a

a

a

a

Ф

11

1s

1r

1m

1

a

s1

a

ss

a

sr

a

sm

Ф

s

A

p

x

p

a

a

a

a

Ф

r1

rs

rr

sm

r

a

a

a

a

m1

ms

mr

mm

Ф

m

a

a

Ф

a

a

11

1m

1

1s

1r

+

Ф

+

Ф

s

r

a

m1

a

mm

Ф

a

a

mr

m

ms

Тогда окончательная система линейных алгебраических уравнений будет сформирована как

a

a

Ф

a

a

11

1m

1

1s

1r

−

Ф

−

Ф

s

r

a

m1

a

Ф

a

a

mm

m

ms

mr

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

Ф

e

(x, y) a

e

+ b

e

x + c

e

через соответствующую базисную функцию

i

i

i

i

параметры a, b и c находятся из системы уравнений

a

ei

+ b

ei

x

+ c

ei

y

e

i

−1

Ф

a

x

i

y

i

i

i

i

i

1

Ф

e

e

e

e

a

i x

+ c

i y

x

y

Ф

i

+ b

b

i

1

j

Ф

.

j

e

e

j

e

j

i

x

j

y

j

c

e

1

Ф

Ф

a

i

+ b

i x

+ c

i y

k

k

k

k

k

k

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

y

. Где

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

38

параллельно уложенных на глубине b стержней длиной l, находящихся друг от друга на расстоянии d (рис. 3). Удельная проводимость земли составляет γ=0,1 См/м.

b

d

d

Рис. 3. Схема заземлителя Задачи расчета стационарного электрического поля аналогичны задачам

электростатического поля, если заменить

−диэлектрик с диэлектрической проницаемостью r 0 на проводящую среду с удельной проводимостью γ,

−электрическую емкость С на проводимость G,

−объемный заряд ρ на плотность тока j, заряд q на ток i.

Полезные соотношения для расчета сопротивления заземлителя:

.

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

1.Для чего предназначено защитное заземление, в чем принцип работы заземлителя?

2.Каким образом конечные элементы вводятся в машину?

3.Что такое базисная функция?

4.Как составить уравнения по методу конечных элементов?

5.Назовите основные этапы расчета электростатического (стационарного электрического) поля методом конечных элементов.

							4. Рабочее задание
1.	Начертите модель заземлителя на плоскости. Определите рабочую область
	(см. лабораторную работу 1). В данной работе она должна быть квадратной.
	Разбейте область на треугольники в соответствии с рис. 1. Шаг
	дискретизации должен быть не более 0,5b. Пронумеруйте точки, начиная с
	нижнего левого угла, отсчитывая по горизонтали. Размеры заземлителя
	приведены в таблицах 1-2, длина заземлителя l=1м.
	Таблица 1
Номер группы		5, 11, 17				6, 12, 18				1, 7, 13					2, 8, 14			3, 9, 15			4, 10, 16
d, м		1				1,2				1,4					1,6			1,8			2
	Таблица 2
Номер бригады				1, 8, 15			2, 9, 16			3, 10, 17			4, 11, 18			5, 12, 19			6, 13, 20			7, 14, 21
(номер компьютера)
b, м				0,8			1			1,2			1,4			1,6			1,8			2
2.	Запишите граничные условия.
3. Напишите				список из				первых			четырех				конечных элементов с						указанием
	вершин i, j, k треугольников. Нумерация треугольников должна
	соответствовать рис. 1.
4. Запишите				выражение				для			расчета				параметров			базисных				функций
	треугольников. Выведите выражения для расчета составляющих векторов
	напряженности электрического поля и плотности тока по оси x и y.
5. В		вычислительной						среде			MATLAB				используйте				шаблон			расчета
	стационарного поля заземлителя методом конечных элементов. Скопируйте
											39

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

%Введите исходные данные

U = ____________;

h = ____________; % шаг дискретизации gamma = ____________;

%Используется сетка (n+1)*(n+1). Укажите ваше n n=____________;

%Через запятую укажите точки с нулевым потенциалом. Программа рассчитана

%на то, что нумерация начинается с точки в левом нижнем углу, отсчет ведется по

%горизонтали.

%Пример: points_0 = [1, 2, 3];

points_0 = [____________];

%Через запятую в первом ряду укажите точки, для которых требуется

%выполнение однородного условия Неймана. Во втором ряду укажите точки,

%потенциал которых должен быть равен потенциалу в точках из первого ряда.

%Пример: points_dfdn_0 = [1, 4, 7;

%

2, 5, 8];

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

%во втором ряду. Тем самым будет выполнено однородное условие Неймана points_dfdn_0 = [____________;

____________];

%Укажите точки с заданным потенциалом U = 1000В. Пример: points_U = [89,93]; points_U = [____________];

number_of_points = (n+1)^2; number_of_elements = 2*n^2; % Определение типа точек

point_type = zeros(number_of_points,1); %Внутренние точки point_type(points_dfdn_0(1,:)) = 2; %Граничные точки с нулевым условием Неймана point_type([points_U, points_0]) = 1; %Граничные точки с заданным потенциалом point_type2_pair = zeros(number_of_points,1);

point_type2_pair(points_dfdn_0(1,:)) = points_dfdn_0(2,:)';

ijk_element = zeros(number_of_elements,3); for q=1:n

ijk_element(n*(q*2-2)+1:n*(2*q-1),1)=(q*(n+1)+1:q*(n+1)+n); ijk_element(n*(q*2-2)+1:n*(2*q-1),2)=(1+(q-1)*(n+1):1+(q-1)*(n+1)+n-1); ijk_element(n*(q*2-2)+1:n*(2*q-1),3)=(q*(n+1)+2:q*(n+1)+n+1);

ijk_element(n*(q*2-2)+n+1:n*(2*q-1)+n,1)=ijk_element(n*(q*2-2)+1:n*(2*q-1),2)+1; ijk_element(n*(q*2-2)+n+1:n*(2*q-1)+n,2)=ijk_element(n*(q*2-2)+1:n*(2*q-1),3); ijk_element(n*(q*2-2)+n+1:n*(2*q-1)+n,3)=ijk_element(n*(q*2-2)+1:n*(2*q-1),2);

end

% Формирование матрицы А

A = sparse(number_of_points,number_of_points);

40

for q=1:number_of_elements i=ijk_element(q,1); j=ijk_element(q,2); k=ijk_element(q,3);

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

%Введите формулы расчета вклада каждого треугольника в матрицу А

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

A(i,i)=____________; A(i,j)=____________; A(j,i)=____________; A(i,k)=____________; A(k,i)=____________; A(j,j)=____________; A(k,k)=____________;

end

%Метод простой итерации

phi_new = zeros(number_of_points,1); phi_new(points_U) = U;

epsilon = 1; itr=0;

while epsilon>0.001 phi_old = phi_new;

for i=1:number_of_points switch point_type(i)

case 0

phi_new(i) = - A(i,[1:i-1 i+1:end])*phi_old([1:i-1 i+1:end])/A(i,i); case 2

phi_new(i) = phi_old(point_type2_pair(i));

end end

epsilon = norm(phi_new-phi_old); itr = itr +1;

if round(itr/100)*100==itr

fprintf('Шаг %d, точность %f\n', itr, epsilon) end

end

phi = phi_new; % потенциалы в точках

Phi = rot90(reshape(phi,n+1,n+1)); % потенциалы в пространстве

% Расчет параметров базисных функций hx = 0:h:n*h;

hy = 0:-h:-n*h;

[X,Y] = meshgrid(hx,hy); %пространственные матрицы координат x и y

y = reshape(rot90(rot90(Y')),number_of_points,1); % координаты y точек рабочей области x = reshape(X',number_of_points,1); % координаты x точек рабочей области

abc_element = zeros(number_of_elements,3); Ex_element = zeros(number_of_elements,1); Ey_element = zeros(number_of_elements,1);

41

Jx_element = zeros(number_of_elements,1); Jy_element = zeros(number_of_elements,1); for q=1:number_of_elements

i=ijk_element(q,1); j=ijk_element(q,2); k=ijk_element(q,3);

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

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

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

abc_element(q,:) = (inv([____________;...

____________;...

____________;])...

* [____________; ...

____________; ...

____________])';

a = abc_element(q,1);

b= abc_element(q,2);

c= abc_element(q,3);

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

%Введите формулы для расчета составляющих векторов напряженности

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

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

Ex_element(q) = ____________; Ey_element(q) = ____________; Jx_element(q) = ____________; Jy_element(q) = ____________;

end

% Построение картины поля

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

Ey = Eyneg; hx = 0:h:n*h; hy = -(0:h:n*h);

contour(hx,hy,Phi), colorbar, hold on, quiver(hx,hy,Ex,Ey), hold off

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

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

b.Введите информацию о сетке конечных элементов и номерах граничных точек;

c.Введите формулы для формирования вклада каждого треугольника в матрицу А;

d.Введите выражение для расчета параметров базисной функции треугольника;

e.Введите формулы для расчета составляющих векторов напряженности электрического поля и плотности тока по оси x и y

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

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

42

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

0					250
-1
-2					200
-3
-4
					150
-5
-6
-7					100
-8
-9	2	4	6	8	50
0

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

7.Запишите вклад в матрицу A первых двух треугольников.

8.Определите потенциал в точке (7·b/3,-15·d/7). Для этого определите, к какому треугольнику относится эта точка, найдите параметры его базисной функции. Сравните результат с расчетом в программе FlexPDE.

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

10.Ток растекается от заземлителя к бесконечно удаленной границе с нулевым потенциалом. Он будет одним и тем же, если его рассчитать через любую замкнутую поверхность, которая охватывает стержни. Это утверждение аналогично теореме Остроградского-Гаусса для электростатического поля. Оно также справедливо согласно 1-му закону Кирхгофа, если стержни заземлителя представить в виде узла с подключенным источником тока.

Для расчета тока необходимо выбрать замкнутую поверхность, охватывающую стержни. Удобно использовать параллелепипед глубиной l (проекция на плоскую область расчета – прямоугольник).

Проведите прямоугольник, охватывающий стержни, и определить нормальные составляющие векторов плотности тока к сторонам

прямоугольника (к горизонтальным сторонам – составляющие Jy, к

вертикальным – Jx). Рассчитайте ток растекания, сложив токи сквозь элементарные поверхности h·l1. Сравните полученный ток с током, рассчитанным при помощи FlexPDE.

11.Определите сопротивление заземлителя.

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

1 Учтите, что ток растекания во всей области задачи больше, чем в расчетной области.

43