fixed width

3.1 Vlnovody

Vývoj programu

Obr. 3.1D.1
Obr. 3.1D.1Příklad dvojrozměrné sítě konečných prvků. Modrá čísla: lokální uzly. Červená čísla: globální uzly.

Metoda konečných prvků je obecná numerická metoda pro řešení parciálních diferenciálních rovnic. Jelikož Maxwellovy rovnice, které popisují šíření elektromagnetické vlny ve vlnovodu, patří do rodiny parciálních diferenciálních rovnic, můžeme metodu konečných prvků k analýze pole ve vlnovodu použít.

Základní postup při výpočtu rozložení pole ve vlnovodu lze shrnout do několika kroků:

  1. Profil vlnovodu rozdělíme na konečné prvky. V případě podélně homogenního vlnovodu rozdělíme profil na malé trojúhelníky. Každý trojúhelník je konečný prvek, vrcholy trojúhelníků nazýváme uzly.
  2. Rozložení podélné složku elektrické intenzity u vidů TM a magnetické intenzity u vidů TE formálně aproximujeme nad každým konečným prvkem lineární, kvadratickou nebo kubickou funkcí. Formální aproximací rozumíme skutečnost, že funkce je vyjádřena pomocí neznámých hodnot intenzity v uzlech a pomocí známých bázových funkcí (lineárních, kvadratických, kubických). Každá bázová funkce nabývá jednotkové hodnoty v jednom uzlu a nulové hodnoty v uzlech ostatních. Aproximace je formálně jednou rovnicí pro N neznámých uzlových hodnot.
  3. Formální aproximaci dosadíme do výchozího vztahu, tj. do řešené rovnice. Jelikož aproximace díky své přibližnosti nesplňuje řešenou rovnici zcela přesně, musíme k našemu vztahu přičíst chybovou funkci (tzv. reziduum). Čím menší je reziduum, tím přesnější je aproximace.
  4. Reziduum minimalizujeme metodou vážených reziduí. Tato metoda spočívá v postupném vynásobením rezidua vhodnými váhovými funkcemi a v integrování tohoto součinu přes celý profil analyzovaného vlnovodu. Zvolíme-li za váhové funkce bázové funkce všech uzlů, v nichž neznáme hodnotu intenzity, dostaneme N rovnic pro N neznámých. Popsaná volba váhových funkcí je nazývána Galerkinovou metodou.
  5. Vyřešíme maticovou rovnici, kterou dostaneme jako výsledek Galerkinovy metody. Tím získáme hodnoty doposud neznámých uzlových intenzit. Maticová rovnice, kterou řešíme v případě vlnovodu, je tzv. zobecněný problém vlastních čísel. Výsledkem je vektor vlastních čísel (vektor kvadrátů kritických vlnových čísel jednotlivých vidů) a matice odpovídajících vlastních vektorů (uzlové hodnoty podélných složek intenzit jednotlivých vidů).
  6. Dosazením uzlových hodnot do formální aproximace získáme aproximační funkci hledané podélné složky intenzity v každém bodě profilu analyzovaného vlnovodu. Z podélné složky intenzity můžeme vypočíst všechny ostatní složky vektorů intenzity elektrického i magnetického pole.

Při praktických výpočtech metodou konečných prvků postupujeme poněkud odlišně. V literatuře jsou totiž k dispozici již hotové matice koeficientů pro normovaný konečný prvek. Stačí nám vzít tyto matice koeficientů, upravit je pro naše konkrétní prvky a spojit tyto matice společně s konečnými prvky do sítě. Celou proceduru lze opět popsat v několika krocích:

  1. Profil vlnovodu rozdělíme na konečné prvky. Tento úkol je naprosto shodný s předchozím postupem.
    Obr. 3.1D.2
    Obr. 3.1D.2Dvojrozměrný konečný prvek pro lineární aproximaci
  2. V literatuře nalezneme vztahy pro výpočet koeficientů jednotlivých konečných prvků. Dosazením plochy e-tého konečného prvku A(e) a úhlů u vrcholů e-tého trojúhelníkového prvku (viz obr. 2) do těchto vztahů dostaneme koeficienty pro naše konkrétní prvky. Vztahy pro matice koeficientů lineární aproximace následují:
    S(e)=n=13Qncotgϑn(e),   T(e)=A(e)12[ 211121112 ],
    Q1=12[ 0000+1101+1 ],   Q2=12[ +10100010+1 ],   Q3=12[ +1101+10000 ]
    ( 3.1D.1 )
  3. Sestavíme diagonální matice pro "izolované" konečné prvky, tj. pro prvky dosud nespojené do sítě (symboly 0 značí nulové matice o rozměrech 3x3)
    S=[ S(1)S(2)S(3)... ],   T=[ T(1)T(2)T(3)... ] ( 3.1D.2 )
  4. Sestavíme vazební matici, která určuje vzájemný vztah mezi lokálními uzly (modrá čísla v obr. 1) a uzly globálními (červená čísla v obr. 1). Sloupce matice odpovídají globálním uzlům, řádky matice uzlům lokálním. Pokud má být lokální uzel přidružen k příslušnému uzlu globálnímu, napíšeme do řádku odpovídajícího číslu lokálního uzlu ve sloupci odpovídajícím uzlu globálnímu hodnotu 1. V opačném případě je hodnota prvku matice nulová. Pro uzly z obr. 1 by vazební matice vypadala následovně:
    Obr. 3.1D.3
    Obr. 3.1D.3Vazební matice
  5. Sdružíme nezávislé konečné prvky do propojené sítě elementů. Matematicky tomu odpovídají následující maticová operace:
    Sc=CTSC,   Tc=CTTC ( 3.1D.3 )
  6. Vyřešením maticové rovnice (je stejná jako u dříve popsaného postupu)
    SE+k2TE=0 ( 3.1D.4 )
    získáme kritická vlnová čísla jednotlivých vidů k a uzlové hodnoty intenzity elektrického pole E těchto vidů. Dosazením uzlových hodnot do formální aproximace získáme aproximační funkci hledané podélné složky intenzity v každém bodě profilu analyzovaného vlnovodu.

Náš matlabovský program je založen na právě popsaném algoritmu. Zjednodušený výpis matlabovského kódu následuje. Jelikož kód je pečlivě komentován, nebudeme ho dále nijak popisovat.

function kn = linearTE( Nx, Ny);

iNx = Nx + 1;
iNy = Ny + 1;

a = 22.86e-3; % šířka vlnovodu
b = 10.16e-3; % výška vlnovodu

dx = ones(1,Nx) * (a/Nx); % vektory rozměrů k.p.
dy = ones(1,Ny) * (b/Ny);

Q1 = [ 0 0 0; 0 1 -1; 0 -1 1] / 2;
Q3 = [ 1 -1 0; -1 1 0; 0 0 0] / 2;

Te = [ 2 1 1; 1 2 1; 1 1 2] /12;

N = 2 * Nx * Ny; % celkový počet k.p.

St = sparse( 3*N, 3*N); % matice pro izolované k.p.
Tt = sparse( 3*N, 3*N);

n = 0;

for ny=1:Ny
for nx=1:Nx
n = n + 1;
lw = 3*n-2;
hg = 3*n;
St(lw:hg,lw:hg) = Q1 * dx(nx)/dy(ny) + Q3 * dy(ny)/dx(nx);
Tt(lw:hg,lw:hg) = Te * dx(nx)*dy(ny)/2;
St(lw+3*Nx:hg+3*Nx,lw+3*Nx:hg+3*Nx) = St(lw:hg,lw:hg);
Tt(lw+3*Nx:hg+3*Nx,lw+3*Nx:hg+3*Nx) = Tt(lw:hg,lw:hg);
end
n = n + Nx;
end

C = get_c1( Nx, Ny, N)

S = C'*St*C; % sdružení izolovaných k.p.
T = C'*Tt*C;

clear C St Tt

[H,K] = eig( full(S), full(T)); % řešení problému vlast.čísel

kn = sqrt( diag( K)); % vektor kritických vlnových čísel

Copyright © 2010 FEEC VUT Brno All rights reserved.