5.1 Modelování drátových antén v časové oblasti metodou momentů

Vývoj programu

Program pro analýzu drátového dipólu implicitním přístupem je z hlediska implementace náročný. Nejdříve je nutné vypočítat budicí Gaussův impuls modulovaný harmonickým signálem (5.1A.7):

EG=E0*4/(c*T*sqrt(pi)).*exp(-(4/T)^2*(t-t0-delt/2).*(t-t0-delt/2)).*cos(2*pi*f0*(t-delt/2));

V druhém kroku vyčíslíme hodnoty koeficientů κ( m, n) a κ( m±, n±) podle vztahu (5.1B.5) resp. (5.1B.13) a dále ještě vypočítáme vzdálenost mezi diskretizačními elementy dle (5.1B.7) a (5.1B.12). Máme možnost volby pro tento výpočet. Buď můžeme použít přesné vztahy (5.1B.7a) resp. (5.1B.12a), nebo přibližně se zanedbáním poloměru anténního vodiče dle (5.1B.7b) a (5.1B.12b) pokud cΔt>>a.

Rmn=zeros(Nz+1,1);
kappa=zeros(Nz+1,1); % priprava vypoctu kappa dle vztahu (5.1B.5)
for m=0:Nz % vztah (5.1B.13)
hlp1=log((m+0.5)*delz + sqrt(((m+0.5)*delz)^2+a^2));
hlp2=log((m-0.5)*delz + sqrt(((m-0.5)*delz)^2+a^2));
kappa(m+1)=(hlp1-hlp2);
% vzdalenost diskret. elem. (5.1B.7) a (5.1B.12), mame moznost volby (5.1B.7a) a (5.1B.12a)
% pokud je splneno c*delt>>a je lepe pouzit (5.1B.7b) a (5.1B.12b) a zanedbat polomer vodice (viz Vrstva B)
% R(m+1)=sqrt(a^2+(m*delz)^2); % (5.1B.7a) a (5.1B.12a)
R(m+1)= m*delz; % (5.1B.7b) a (5.1B.12b)
end
kappaFI=1/4/pi/ep*kappa; % kappa pro skalarni pot. (5.1B.13)
kappa=mi/4/pi*kappa; % kappa pro vektorovy pot. (5.1B.5)

Dále vyčíslíme hodnoty prvků matice vystupujících v rovnici (5.1B.35), tj. [A11C(m, n)], [A12C(m, n)] dle (5.1B.35) a (5.1B.22) pro vektorový potenciál a [φ11C(m, n)], [φ12C(m, n)] a [φ13C(m, n)] dle (5.1B.36), (5.1B.27b) a (5.1B.29) pro skalární potenciál. Zde je nutné znovu upozornit, že členy těchto matic se týkají proudů v čase tRk>tk-1. V případě matic vektorového potenciálu [A11C(m, n)], [A12C(m, n)] je zdrojový text následující:

A11C=zeros(N,N); % koef. vekt. pot. dle (5.1B.35) a (5.1B.22)
A12C=zeros(N,N);
for m=1:N,
for n=1:N,
mn=abs(m-n); % uvazuj jen cleny v case tk>=tRk>tk-1
if(R(mn+1)<cdelt) % tj. vzdalenst musi byt mensi nez c*delt
A11C(m,n)=A11C(m,n)+(1-R(mn+1)/cdelt)*kappa(mn+1);
A12C(m,n)=A12C(m,n)+R(mn+1)/cdelt*kappa(mn+1);
end
end
end

V případě matic skalárního potenciálu [φ11C(m, n)], [φ12C(m, n)] a [φ13C(m, n)] uvedeme zde jen část kódu pro výpočet příspevku prvního členu s indexy (m+, n+) v (5.1B.36), výpočet ostatních členů s indexy (m+, n-), (m-, n+) a (m-, n-) je analogický:

FI11C=zeros(N,N); % koef. skal. pot. dle (5.1B.36)
FI12C=zeros(N,N); % pro ++ dle (5.1B.27b) a (5.1B.29);
FI13C=zeros(N,N); % pro ostatni kombinace analogicky

for m=1:N,
for n=1:N,
mn=abs(m-n); % vypocet indexu vzdalenosti mezi elementy
mnPP=mn; % P je +; M je -
mnPM=mn+1;
mnMP=abs(mn-1);
mnMM=mn;
if(m<n)
hlp=mnPM;
mnPM=mnMP;
mnMP=hlp;
end
% prispivaji jen cleny, pro ktere plati tRk>tk-1
if(R(mnPP+1)<cdelt)
FI11C(m,n)=FI11C(m,n) -delt/delz/2*(1-R(mnPP+1)/cdelt)^2*kappaFI(mnPP+1); %(5.1B.29a)
FI12C(m,n)=FI12C(m,n) -delt/delz/2*(1-(R(mnPP+1)/cdelt)^2)*kappaFI(mnPP+1); %(5.1B.29b)
FI13C(m,n)=FI13C(m,n)-1/delz*kappaFI(mnPP+1); %(5.1B.27b)
end
.
.
.
end
end

Protože inverzní matice v soustavě (5.1B.35) nezávisí na čase, vyčíslíme ji před časovou smyčkou:

invA11CFI11C=inv(A11C+delt/delz/2*FI11C); % vypocet inverzni matice dle (5.1B.35) leva strana

Nakonec implementujeme časovou smyčku, v níž počítáme na základě známých časových vzorků vektorového potenciálu a skalárního potenciálu vzorky budoucí. Poznamenejme, že hodnoty proudů, nebo hodnoty integrálů z proudů dle času, které nepadnou do diskrétních okamžiků, jsou mezi těmito okamžiky lineárně aproximované. V proměnné IS je uložen integrál proudu, který je počítán dle lichoběžníkového pravidla. Význam proměnných by měl být zřejmý. Opět zde v případě výpočtu skalárního potenciálu uvádíme zdrojový kód jen pro výpočet příspěvku prvního členu s indexy (m+, n+) viz (5.1B.34) a (5.1B.36), vypočet ostatních příspěvků členů s indexy (m+, n-), (m-, n+) a (m-, n-) je analogický.

% CASOVA SMYCKA
for k=3:Nt,
A2=zeros(N,1); % (5.1B.19b)
FI2=zeros(N,1); % (5.1B.36d), (5.1B.25b) a dale analog. pro +-,-+,--
for m=1:N,
for n=1:N,
% vypocet indexu vzdalenosti mezi elementy
mn=abs(m-n); % vektorovy potencial
mnPP=mn; % skalarni potencial
mnPM=mn+1; % P je +; M je -
mnMP=abs(mn-1);
mnMM=mn;
if(m<n)
hlp=mnPM;
mnPM=mnMP;
mnMP=hlp;
end
% dale jsou pocitany jen cleny, pro ktere plati: tRk<=tk-1
if(R(mn+1)>=cdelt) % vektorovy potencial (5.1B.19b)
tk=k-R(mn+1)/cdelt;
if(tk>1)
ftk=floor(tk); % linearni aproximace proudu
In=(I(n,ftk+1)-I(n,ftk))*(tk-ftk)+I(n,ftk);
A2(m,1)=A2(m,1)+In*kappa(mn+1);
end
end
if(R(mnPP+1)>=cdelt) % skalarni potencial ++; (5.1B.36d), (5.1B.25b)
tk=k-R(mnPP+1)/cdelt;
if(tk>1)
ftk=floor(tk); % linearni aproximace integralu proudu
ISn=(IS(n,ftk+1)-IS(n,ftk))*(tk-ftk)+IS(n,ftk);
FI2(m,1)=FI2(m,1)+ISn/delz*kappaFI(mnPP+1);
end
end
.
.
.
end
end
E(nap,1)=EG(k); % vektor buzení je nenulovy jen v miste napajeni (5.1B.7)
% (5.1B.35)
I(:,k)=invA11CFI11C*(delt*E-A12C*I(:,k-1)-A2+A -delt/2/delz*(FI12C*I(:,k-1)+FI13C*IS(:,k-1)+FI2+FI)); % integrace proudu dle lichobeznikoveho pravidla (5.1B.15)
IS(:,k)=(I(:,k)+I(:,k-1))*delt/2+IS(:,k-1); % vypocet pro dalsi krok
A=A11C*I(:,k)+A12C*I(:,k-1)+A2; % (5.1B.18) a (5.1B.21) FI=FI11C*I(:,k)+FI12C*I(:,k-1)+FI13C*IS(:,k-1)+FI2; % (5.1B.30)-(5.1B.33)
% s uvazenim (5.1B.26)-(5.1B.29)
end

Copyright © 2010 FEEC VUT Brno All rights reserved.