MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Математика\Partial Differential Equations (PDE) Toolbox

В.Е.Шмелев "Примеры применения 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 (пФ/мм)

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


Поиск по сайту:

Система Orphus

Яндекс.Метрика