В.Е.Шмелев "Примеры применения PDE Toolbox"

Расчёт электростатического поля “коаксиального” кабеля со смещённой жилой

Пусть имеется цилиндрическая электролитическая ванна, моделирующая “коаксиальный” кабель со смещённой жилой. Пусть радиус оболочки rob=250 мм, радиус жилы rz=20 мм, смещение осей оболочки и жилы d=150 мм, напряжение U=10 В, относительная диэлектрическая проницаемость epsilon=1. Расчёт эквипотенциалей проведём аналитически и в PDE Toolbox. Ниже приведен текст вычислительного сценария аналитического расчёта.

% vannaks - Расчёт электростатического полq в "коаксиальном" кабеле со смещённой жилой
%
% Входные данные: epsilon - проницаемость;
% rob - радиус оболочки; rz - радиус жилы;
% d - смещение оси жилы относительно оси оболочки;
% U - напрqжение;
% nf - число шагов по потенциалу.
%
% Выходные данные: c0 - ёмкость на единицу длины;
% rk - радиусы эквипотенциалей.
%
% В обычной фигуре строитсq картина эквипотенциалей
%
eps0=8.85419e-3; % Абсолютнаq диэлектрическаq проницаемость вакуума, пФ/мм
if exist('epsilon','var'), sepsilon=num2str(epsilon); else sepsilon='1'; end
if exist('rob','var'), srob=num2str(rob); else srob='250'; end
if exist('rz','var'), srz=num2str(rz); else srz='20'; end
if exist('d','var'), sd=num2str(d); else sd='40'; end
if exist('U','var'), sU=num2str(U); else sU='10'; end
if exist('nf','var'), snf=num2str(nf); else snf='10'; end
SS=inputdlg({'epsilon','rob (миллиметры)','rz (миллиметры)','d (миллиметры)','U (вольты)','nf шагов по потенциалу'},...
'Ввод исходных данных',1,{sepsilon,srob,srz,sd,sU,snf});
epsilon=eval(SS{1}); rob=eval(SS{2}); rz=eval(SS{3}); d=eval(SS{4}); U=eval(SS{5}); nf=eval(SS{6});
disp(['epsilon=',num2str(epsilon),'; rob=',num2str(rob),'; rz=',num2str(rz),'; d=',num2str(d),'; U=',num2str(U),'; nf=',num2str(nf)])
s1=(rob^2-rz^2-d^2)/2/d;
s2=(rob^2-rz^2+d^2)/2/d;
a=sqrt(s1^2-rz^2);
c0=2*pi*eps0*epsilon/log((s2-a)*(s1+a)/rob/rz)
tau=c0*U;
fz=tau*log((s1+a)/rz)/(2*pi*eps0*epsilon);
fob=tau*log(rob/(s2-a))/(2*pi*eps0*epsilon);
fk=linspace(0,U,nf+1);
hi=((s2-a)*(s1+a)/rob/rz).^((fob+fk)/U);
x=s2-a*(hi.^2+1)./(hi.^2-1)
rk=2*a*abs(hi./(1-hi.^2))
t=0:0.004*pi:2*pi;
for k=1:nf+1
    plot(rk(k)*cos(t)+x(k),rk(k)*sin(t),'k-')
    hold on
end
grid on

В PDE Toolbox рассчитаем не только распределение потенциала, но и линейную плотность заряда жилы, энергию электрического поля, ёмкость кабеля на единицу длины. Последний интегральный параметр будем определять двумя способами: 1) через двумерный интеграл (через энергию) и 2) через криволинейный интеграл (через поток вектора электрического смещения на единицу длины кабеля), причём интеграл вычислим для нескольких контрольных замкнутых кривых.

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

Рис.1. Поперечное сечение кабеля (ванны)

На рис.2 показана картина аналитически рассчитанных эквипотенциалей.

Рис.2. Картина аналитически рассчитанных эквипотенциалей

Выходная информация в командном окне:

>>vannaks
epsilon=1; rob=250; rz=20; d=150; U=10; nf=10
c0 =
    0.026828 (пФ/мм)
x =
  Columns 1 through 6
    5.6843e-014    67.921    101.72    120.63    131.91    138.89
  Columns 7 through 11
    143.31    146.16    148    149.21    150 (мм)
rk =
  Columns 1 through 6
    250    169.72    124.41    94.954    74.189    58.794
  Columns 7 through 11
    47.011    37.803    30.51    24.683    20 (мм)

Теперь прорисуем геометрию, изображённую на рис.1, в GUI-приложении PDETool:

pdecirc(0,0,250,'Cob') % Оболочка
pdecirc(150,0,20,'Cz') % Жила
pdecirc(x(2),0,rk(2),'C1') % Эквипотенциаль 1 В
pdecirc(x(3),0,rk(3),'C2') % Эквипотенциаль 2 В
pdecirc(x(4),0,rk(4),'C3') % Эквипотенциаль 3 В
pdecirc(x(5),0,rk(5),'C4') % Эквипотенциаль 4 В
pdecirc(0,0,185,'CC1')
pdecirc(0,0,220,'CC2')

Потенциал оболочки зададим равным нулю, потенциал жилы 10 В. Граничные окружности геометрических объектов C1, C2, C3, C4, CC1, CC2 используем в качестве линий интегрирования. Технологию задания граничных условий, коэффициентов PDE, генерации сетки и запуска решения PDE показывать не будем. На рис.3 показана разметка зон расчётной области. На рис.4 показана картина эквипотенциальных линий, рассчитанная в PDETool.

Рис.3. Зоны расчётной области

Рис.4. Картина эквипотенциальных линий, рассчитанная в PDETool

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

[ux,uy]=pdegrad(p,t,u);
ar=pdetrg(p,t);
energy=eps0*epsilon*(ux.^2+uy.^2)*ar.’/2
energy =
    1.341 (пДж/мм)
c0p=2*energy/U^2
c0p =
    0.026819 (пФ/мм)

Как видно, значения c0 и c0p совпадают с точностью до трёх верных значащих цифр. Вычисления проводились на сетке с числом узлов 3343 и числом элементов 6532.

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

; ; ;

.

Приведённые соотношения записаны для случая, когда положительное направление интегрирования – против часовой стрелки, и положительная нормаль – внешняя. Теперь можно кодировать последовательность операторов MATLAB:

Dx=-ux*eps0*epsilon;
Dy=-uy*eps0*epsilon;
Dxu=pdeprtni(p,t,Dx);
Dyu=pdeprtni(p,t,Dy);
% Интегрируем по границе жилы.
% Граничные элементы, у которых зона 8 справа,
% будем учитывать с неизменным знаком, остальные – с противоположным.
iep=find(e(7,:)==8&e(6,:)==0) % индексы положительных граничных элементов
iem=find(e(6,:)==8&e(7,:)==0) % индексы отрицательных граничных элементов
Dxp=(Dxu(e(1,iep))+Dxu(e(2,iep)))/2;
Dxm=(Dxu(e(1,iem))+Dxu(e(2,iem)))/2;
Dyp=(Dyu(e(1,iep))+Dyu(e(2,iep)))/2;
Dym=(Dyu(e(1,iem))+Dyu(e(2,iem)))/2;
dxp=p(1,e(2,iep))-p(1,e(1,iep));
dxm=p(1,e(2,iem))-p(1,e(1,iem));
dyp=p(2,e(2,iep))-p(2,e(1,iep));
dym=p(2,e(2,iem))-p(2,e(1,iem));
tauz=dyp*Dxp-dxp*Dyp-dym*Dxm+dxm*Dym % Интеграл по границе жилы готов

Чтобы вычислить аналогичные интегралы для других кривых, последовательность операторов от iep=… до tauz=… нужно выполнить много раз для различных пар зон. Значит, эту последовательность целесообразно оформить в виде функции MATLAB с входными параметрами p, e, Dxu, Dyu, zonel, zoner и с одним выходным параметром tau:

function tau=intpotok(p,e,Dxu,Dyu,zonel,zoner);
iep=find(e(7,:)==zoner&e(6,:)==zonel); % индексы положительных граничных элементов
iem=find(e(6,:)==zoner&e(7,:)==zonel); % индексы отрицательных граничных элементов
Dxp=(Dxu(e(1,iep))+Dxu(e(2,iep)))/2;
Dxm=(Dxu(e(1,iem))+Dxu(e(2,iem)))/2;
Dyp=(Dyu(e(1,iep))+Dyu(e(2,iep)))/2;
Dym=(Dyu(e(1,iem))+Dyu(e(2,iem)))/2;
dxp=p(1,e(2,iep))-p(1,e(1,iep));
dxm=p(1,e(2,iem))-p(1,e(1,iem));
dyp=p(2,e(2,iep))-p(2,e(1,iep));
dym=p(2,e(2,iem))-p(2,e(1,iem));
tau=dyp*Dxp-dxp*Dyp-dym*Dxm+dxm*Dym; % Интеграл по границе раздела двух зон

Итак, интегрируем по всем кривым:

Dx=-ux*eps0*epsilon;
Dy=-uy*eps0*epsilon;
Dxu=pdeprtni(p,t,Dx);
Dyu=pdeprtni(p,t,Dy);
% Граница жилы:
tauz=intpotok(p,e,Dxu,Dyu,0,8)
% Эквипотенциаль 4 В:
tau4v=intpotok(p,e,Dxu,Dyu,8,6)+intpotok(p,e,Dxu,Dyu,13,12)
% Эквипотенциаль 3 В:
tau3v=intpotok(p,e,Dxu,Dyu,6,2)+intpotok(p,e,Dxu,Dyu,12,11)
% Эквипотенциаль 2 В:
tau2v=intpotok(p,e,Dxu,Dyu,2,1)+intpotok(p,e,Dxu,Dyu,11,10)+intpotok(p,e,Dxu,Dyu,11,7) +intpotok(p,e,Dxu,Dyu,14,9)
% Эквипотенциаль 1 В:
tau1v=intpotok(p,e,Dxu,Dyu,1,3)+intpotok(p,e,Dxu,Dyu,10,4)+intpotok(p,e,Dxu,Dyu,7,4) +intpotok(p,e,Dxu,Dyu,9,5)
% Граничнаq окружность CC1:
tau1=intpotok(p,e,Dxu,Dyu,3,4)+intpotok(p,e,Dxu,Dyu,1,7)+intpotok(p,e,Dxu,Dyu,2,11) +intpotok(p,e,Dxu,Dyu,6,12)+intpotok(p,e,Dxu,Dyu,8,13)+intpotok(p,e,Dxu,Dyu,1,10)
% Граничнаq окружность CC2:
tau2=intpotok(p,e,Dxu,Dyu,4,5)+intpotok(p,e,Dxu,Dyu,7,9)+intpotok(p,e,Dxu,Dyu,9,5) +intpotok(p,e,Dxu,Dyu,10,9)
% Вычисляqм ёмкости.
c0int=[tauz tau4v tau3v tau2v tau1v tau1 tau2]/U

Выходные данные в командном окне:

>> tauz =
    0.22344
tau4v =
    0.26731
tau3v =
    0.26873
tau2v =
    0.26798
tau1v =
    0.26824
tau1 =
    0.2689
tau2 =
    0.34075 (пКл/мм)
c0int =
  Columns 1 through 6
    0.022344 0.026731 0.026873 0.026798 0.026824 0.02689
  Column 7
    0.034075 (пФ/мм)

Из этого примера отчётливо видно, что интегрирование по границам обладает большей численной неустойчивостью, чем интегрирование по двумерной расчётной области или по подобластям.