Svobodný software v rukách vědy

22.9.2008 17:29 | Přečteno: 2588× | Něco off-topic | | poslední úprava: 29.9.2008 08:47

Původně jsem zamýšlel tento zápisek věnovat pouze krátkému skriptu, s jehož pomocí se mi podařilo v bashi jednoduše ořezat obrázky a vytvořit z nich animovaný gif. Jak ale tento úmysl dostával v mé mysli jasnější a jasnější obrysy bylo mi jasné, že by tomuto zápisku jaksi chyběla hlava i pata - bylo by v něm pouze co se povedlo, ale ne proč a jaký to bude mít další efekt. Proto jsem se rozhodl tento zápisek rozšířit až na současnou formu.

---

• Úvod
• FEMM a ti druzí
• Konkrétní aplikace
• Příprava modelu
• LuaScript
• Zpracování číselných dat
• Zpracování obrazového výstupu
• Závěr

---

Úvod

Před rokem jsem byl přijat na doktorské studium na Katedře elektromechaniky a výkonové elektroniky na FEL ZČU v Plzni, kde tou dobou začínal výzkum zabývající se lámáním tyčí rotorového vinutí asynchronních strojů (grant MPO č. 2A-2TP1/139), do kterého jsem nakonec také byl zapojen. Kromě logicky nezbytné části výzkumu prováděné v proprietárním software (ANSYS, MATLAB) jsem jisté části své práce prováděl i v opensource software, zejména v programu FEMM (AFPL).

FEMM a ti druzí

Právě v programu FEMM byly prováděny počáteční výpočty sledovaného asynchronního stroje - důvod je jednoduchý: rychlý preprocessing (usnadněný například možností importu *.dxf výkresů vytvořených v opensource edici programu QCAD (GPL)) a postprocessing a vynikající přesnost, umožňující rychlé získání prvních orientačních výsledků. Do třetice všeho dobrého je možnost skriptování operací v programu FEMM pomocí jazyku Lua (vlastní licence OSS kompatibilní, od verze 5.0 MIT licence). Aby toho nebylo málo, samozřejmostí je zpracování výsledků v OpenOffice.org (LGPL) (přes výstupní soubor) či jejich analýza za pomoci GNU/Octave (GNU/GPL).

Konkrétní aplikace

Protože opensource software již nějaký čas používám a za mnohé mu vděčím, cítím i jistou povinnost mu něco vrátit. Jednou z úloh v rámci programu MPO 2A-2TP1/139 vyřešených čistě za pomoci opensource programů je například studium drážkové harmonické asynchronního stroje. Vzhledem k povaze výsledků jsem zavázán neuvádět konkrétní čísla, nicméně celý postup zveřejnit lze a lze tím i demonstrovat sílu spojeného opensource software s malým bonusem na konci. Uvažujme tedy asynchronní stroj, u nějž se snažíme získat informace o kmitech momentu, které vznikají drážkováním statoru a rotoru. Tyto kmity vznikají jako důsledek vzájemných superpozic drážek stran statoru a rotoru a lze je sledovat pouze tak, že se vždy vypočte elektromagnetické pole v modelu asynchronního stroje pro určitou vzájemnou pozici těchto dvou těles. Zpravidla stačí rotor postupně pootočit o jednu drážkovou rozteč rotoru v rozumném počtu kroků. Co tedy potřebujeme provést abychom se dobrali výsledku:

• vytvořit model v programu FEMM
• vypočítat jej
• pootočit rotorem o určitý úhel
• přepočítat model pro novou polohu

Toto se opakuje v rozumné míře v rozsahu jedné drážkové rozteče.

Příprava modelu
Přípravu modelu v programu FEMM lze rozdělit na dvě části - první, vlastní nákres geometrie (spíše věc software QCAD), a druhou, vlastní definici problému v programu FEMM. To samo o sobě vyžaduje jistou znalost této aplikace, proto se vyhnu detailům a omezím se na konstatování, že v této fázi je nejdůležitějším úkonem přiřazení všech křivek rotoru a definici jednotlivých materiálů jednotné skupině, v tomto konkrétním případě skupině označené číslem "1".

LuaScript
Aby bylo možno sledovat hledanou drážkovou harmonickou, je třeba nějakým způsobem automatizovat dříve uvedené kroky. Dejme tomu, že mezi drážkami rotoru stroje je 10 úhlových stupňů (odpovídá stroji se 36ti drážkami) a tuto vzdálenost chceme projet tak, abychom na ní získali 100 vypočtených modelů. Od našeho LuaScriptu tedy očekáváme, že udělá následující:

• načte vytvořený model
• vybere rotor (skupina "1") a pootočí jej
• uloží upravený model
• vypočte jej
• otevře výsledky v postprocesoru
• označí křivku, podél které se odečítá moment
• odečte moment a vypíše jeho hodnotu do souboru
• vykreslí rozložení magnetického pole ve stroji a obrázek uloží

a to vše zopakuje 100x. Skript schopný provést uvedené úkony pak vypadá následovně:

---

--- otaceni rotorem pro FEMM
--- skript verze 1.0

open("vzor.FEM") --- otevreni vstupniho souboru

vystup=openfile("moment.txt", "w") --- otevreni vystupniho souboru

mi_seteditmode("group") --- v preprocesoru budou editovany skupiny

dAlfa=10/100 --- uhel pootoceni rotoru mezi jednotlivymi kroky

n=0

--- vypocet vzorovaho modelu
mi_saveas(format("motor_%1$d.FEM",n)) --- ulozeni upraveneho modelu
mi_analyse() --- spusteni analyzy

mi_loadsolution() --- nacteni vysledku

mo_showdensityplot(0,0,1,0,"mag") --- zobrazeni indukce ve stroji
mo_savebitmap(format("motor_%1$d.bmp",n)) --- ulozeni vysledku do bitmapy

mo_seteditmode("contour") --- v postprocesoru budou editovany linie
mo_selectpoint(-100,0.5) --- vyber prvniho z bodu na kruznici (horni cast)
mo_selectpoint(100,0.5) --- vyber druheho z bodu na kruznici (horni cast)
mo_selectpoint(-100,-0.5) --- vyber prvniho z bodu na kruznici (dolni cast)
mo_selectpoint(100,-0.5) --- vyber druheho z bodu na kruznici (dolni cast)
MDC,DM,TV,FV=mo_lineintegral(4) --- vypocet momentu stroje na zaklade integralu po linii

write(vystup, MDC, " ", DM, " ", TV, " ", FV, "\n") --- ulozeni vystupnich hodnot

mo_close() --- uzavreni postprocesoru

for n=1,99 do

mi_selectgroup(1) --- vyber skupiny 1 (rotor)
mi_move_rotate(0,0,dAlfa) --- pootoceni rotoru
mi_saveas(format("motor_%1$d.FEM",n)) --- ulozeni upraveneho modelu

mi_analyse() --- spusteni analyzy
mi_loadsolution() --- nacteni vysledku

mo_showdensityplot(0,0,1,0,"mag") --- zobrazeni indukce ve stroji
mo_savebitmap(format("motor_%1$d.bmp",n)) --- ulozeni vysledku do bitmapy

mo_seteditmode("contour") --- v postprocesoru budou editovany linie
mo_selectpoint(-100,0.5) --- vyber prvniho z bodu na kruznici (horni cast)
mo_selectpoint(100,0.5) --- vyber druheho z bodu na kruznici (horni cast)
mo_selectpoint(-100,-0.5) --- vyber prvniho z bodu na kruznici (dolni cast)
mo_selectpoint(100,-0.5) --- vyber druheho z bodu na kruznici (dolni cast)
MDC,DM,TV,FV =mo_lineintegral(4) --- vypocet momentu stroje na zaklade integralu po linii

write(vystup, MDC, " ", DM, " ", TV, " ", FV, "\n") --- ulozeni vystupnich hodnot

mo_close() --- uzavreni postprocesoru

end

closefile(vystup) --- uzavreni vystupniho souboru


---

Pozn.: čísla v LuaScriptu neodpovídají žádnému skutečnému stroji.

Z FEMMu tedy za pomoci uvedeného skriptu získáme dvě skupiny výsledků:
1. výstupní soubor se závislostí velikosti momentu stroje na úhlu mezi statorem a rotorem
2. sérii obrázků zobrazujících magnetické pole ve stroji
Odtud vychází následující zpracování dat.

Zpracování číselných dat
Nejjednodušší možností vizualizace výsledných dat je prosté natažení výsledků do OpenOffice.org, nicméně to je "pouhé" zobrazení. Na skutečnou analýzu dat je třeba použít silnější nástroj - GNU/Octave. Vypočtené hodnoty totiž kromě určitého rozptylu způsobeného metodou řešení obsahují i určitou harmonickou složku, která odpovídá hledané drážkové harmonické. Nezbývá nám tedy než napsat skript pro GNU/Octave, který provede rychlou Fourierovu analýzu vypočtených dat. Tento skript je pak:

---

% Fourierova analyza vypoctenych vysledku
% version 1.0

M_str=1000; % strední hodnota vypocteneho momentu

moment=[ --- zde je vypoctenych 100 hodnot oddelenych carkou --- ]-M_str; % vektor vstupnich dat

uhel=(360/36)/100; % uhel mezi drazkami rotoru

rozvoj=2/length(moment)*fft(moment); % fourieruv rozvoj vektoru vysledku

x=1/uhel*[0:round(length(rozvoj)/2)]/round(length(rozvoj)); % vypocet dat pro osu x

plot(x(1:length(moment)/2),abs(rozvoj(1:round(length(rozvoj)/2)))) % vykresleni vysledku do grafu
title ('Fourierova analýza momentu stroje')
xlabel ('1/alpha')
ylabel ('M [N.m]')

% vypis vysledku do konzole
x(1:length(moment)/2) % osa x (uhlova perioda kmitu)
abs(rozvoj(1:round(length(rozvoj)/2))) % velikost kmitu
arg(rozvoj(1:round(length(rozvoj)/2))) % faze kmitu


---

Po provedení této analýzy již lze jednoduše vybrat maximum (maxima) Fourierova rozvoje momentu a tedy určit drážkové harmonické stroje.

Zpracování obrazového výstupu
Pro účely následných prezentací je třeba skupinu výstupních obrázků nějak upravit a vytvořit z nich (pochopitelně efektní, ale nic skutečně neříkající) animaci. V tom nám ale brání dva zásadní problémy:
1. FEMM ukládá obrázky do *.bmp, který nelze animovat
2. výsledné obrázky jsou obdélníkové a zarovnané doprostřed stránky nevypadají dobře
Rovnou stovku výsledných obrázků je tedy třeba ořezat, konvertovat do animovatelného formátu (gif) a animovat. K tomu nám pomůže již ryze linuxový skript, využívající jednak ImageMagick (licence GPL kompatibilní) (konkrétně utilitu convert) a prográmek gifsicle (GNU/GPL v2). Požadavky na tento skript jsou pak jednoduché - po skriptu, který spustíme v adresáři s obrázky chceme, aby provedl:

• ořezal obrázek na čtverec
• převedl jej do formátu gif
• všechny gify přesunul do jednoho adresáře
• a provedl jejich animaci

V případě výstupních souborů z programu FEMM jde konkrétně na ořezání na rozměr 761x761px a ve formě skriptu vypadá uvedený postup následovně:

---

#!/bin/sh
for f in *.bmp
do convert $f -extent 761x761 $f.gif
done

mkdir gifs

mv *.gif ./gifs/

cd ./gifs

gifsicle -O --delay=10 --loop *.gif > animace.gif


---

Výsledek pak vypadá následovně:

(lze rozkliknout, 13MB)

Závěr

Pomyslným vyvrcholením všeho je opětovný návrat do OpenOffice.org, ve kterých je z uvedeného výzkumu zpracována výzkumná zpráva. Z předloženého příkladu je tedy zřejmé, že za pomoci opensource software lze provést výzkum i tak náročných partií (zde konkrétně z oblasti elektrických strojů) jako je drážková harmonická asynchronního stroje. Veškeré použité prostředky jsou přitom ryze otevřené povahy a v uvedeném textu zdaleka není plně využito jejich možností. Z porovnání s výsledky z uzavřených programů (ANSYS, MATLAB) vyplývá, že při použití uvedeného software nedochází k žádným významným odchylkám ve výsledcích a veškeré výsledky získané za použití uvedených metod jsou důvěryhodné. Je tedy pouze na volbě každého, zda použije uzavřený software nebo se rozhodne využít některou opensource alternativu.

       

Komentáře

[Käyttäjä 11133]

Ona neexistuje možnost jak to popsal analyticky? Nebylo by pro harmonickou analýzu lepší použít software typu Maple/Maxima a udělat to analyticky přes integrály, než s pomocí octave přes sumy? Možná mi něco uniká a ta konstrukce se analyticky popsat nedá, pak se omlouvám, blogpost jsem jen proletěl.

[JiK]

tohle je presne ta uvaha, kde se prozradi, ze jsi fyzik a ne inzenyr.

[Käyttäjä 11133]

No... Řeknu to takhle.. Nevím co jsem, ani nevím co budu... jedno je však jisté... Inženýr to nebude. :D

[pushkin]

Ono to jde i nejde Co lze předpovědět jsou frekvence, na kterých by se ta drážková harmonická měla objevovat, pro tohle vztahy existují. Co už dost dobře předpovědět nejde je amplituda, tu ovlivňuje tvar zubů, který už se analyticky dost dobře popsat nedá, a relace mezi počty zubů statoru a rotoru. Jinak není za co se omlouvat, též dávám přednost ručním výpočtům, pokud je lze provést.

[David Jaša]

Co už dost dobře předpovědět nejde je amplituda, tu ovlivňuje tvar zubů, který už se analyticky dost dobře popsat nedá, a relace mezi počty zubů statoru a rotoru.

Nějaký takový střet s realitou jsem coby student zcela jiného technického oboru čekal.

[Aleš Kapica]

Ad zpracování obrázků, podle mne je mnohem lepší než imagemagick, použít nástrojů z balíku netpbm.

[Daniel Kvasnička ml.]

ZCU? To je to doupe, kde vznikla ta obludna hruznost, ktere rikate IS/STAG?...

[Ladicek]

No počkej, teď tam vznikají ještě lepší věci, klíčové slovo: portál

[pushkin]

Za to můžou FAVácí, my neee

[Ladicek]

To není pravda, CIVácí!

[Daniel Kvasnička ml.]

To je fuk, nejaky obetni beranek se musi najit

[CEST]

Chlape, nemam vubec paru, o cem to pises (jak by tady v okoli rekli "Ich verstehe Bahnhof"), ale pripada mi to jako dobre zpracovana prikladova studie. Co jsem to cetl, tak mi do oka nepadla zadna silena pravopisna chyba a taky je to celkove dobre zpracovany (logicky). Mozna by stalo za to tohle nekam ulozit (at jako clanek nebo jako nejakej dokument treba do root.cz knihovny) pro pripadny zajemce linuxu ve skolstvi nebo naopak jako dukaz pro pripadny odpurce linuxu ve skolstvi, mozna i prumyslu nebo vede.

[otasomil]

Tez pridavam kladne hodnoceni a spickove zpracovani daneho tematu.
Lamani rotorovych tyci v asynchronnich motorech kde se ona klec vytvari litim ve zvlastni kokile pod tlakem je zpusobena naslednym pnutim pri chladnuti a proto se rotorove plechy lisuji nepatrne do "talire" a skladaji stridave proti sobe nacez pri liti klece v kokile nejsou zcela stlaceny (chybi k tomu cca 1 mm u motoru vykonu radove nekolik kW).
Timto se temer zcela odstrani pnuti pri chladnuti rotoru a dojde k optimalnimu stazeni svazku rotorovych plechu.
Toliko praxe povalecnych let kdy se zacala medena svarovana klec nahrazovat hlinikovou strikanou.
Takto to muzete komentovat jako ja praktik v elektrotechnice a elektronice.
Aplikace open source na linuxove platforme velmi potesi obzvlaste ve vede a vyzkumu - jen tak dale.