Файл: Лабораторные работы по ИТ.pdf

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

b

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

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

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

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

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

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

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

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

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

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

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

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

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

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

U = ____________;

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

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

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

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

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

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

points_0 = [____________];

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

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

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

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

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

%во втором ряду. Тем самым будет выполнено однородное условие Неймана 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.

Рис. 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