FlowPro

CFD software written in JAVA


LabyrinthSeal3D is a set of subroutines designated for studying of fluid flow in labyrinth seal. For download click here.

Úvod

Softwarový modul "Labyrinth Seal 3D FlowPro Module" byl vytvořen jako doplňkový výpočtový modul pro otevřený CFD software FlowPro, který je dlouhodobě vyvíjen pracovníky Katedry mechaniky a výzkumného programu 3 Centra excelence NTIS - Nové technologie pro informační společnost na Fakultě aplikovaných věd ZČU v Plzni. Software FlowPro obsahuje výpočetní jádro založené na nespojité Galerkinově metodě konečných prvků a zároveň rozhraní, v němž je možné implementovat vlastní výpočtový modul pro řešení navrženého matematického modelu. Vytvořený výpočtový modul "Labyrinth Seal 3D FlowPro Module" umožňuje provádět numerické simulace proudění páry ve 3D modelech labyrintových rotorových ucpávek parních turbín pro předepsané hodnoty excentricity a úhlové rychlosti precese rotoru. Výpočtový modul obsahuje rovněž programové nástroje pro stanovení silových účinků proudící páry na rotor a pro stanovení frekvenčně závislých matic tuhosti a tlumení ucpávky.

Generování výpočetní síťě

Vytvoření výpočtové síťě je klíčovou činností prováděnou před začátkem výpočtu. Vytvoření kvalitní výpočtové síťě je však časově velmi náročné. Z tohoto důvodu byly vytvořeny skripty v Matlabu pro automatické generování výpočetní sítě. Jako nejvýhodnější se jeví implementace blokově strukturované 2D sítě, která je následně orotována okolo x-ové osy a tím je vytvořena 3D síť. Jednotlivé bloky jsou definovány ve skriptu blockDef.m viz příklad níže. Po spuštění funkce showPoints.m se jednotlivé bloky zobrazí na obrazovce a je možné provést kontrolu sítě, viz obrázek 1.

function [p,TP,N,z] = blockDef(e) % vrcholy p = []; p = [p;[0,0]]; %1 p = [p;[1,0]]; %2 p = [p;[1,1]]; %3 p = [p;[0,1]]; %4 p = [p;[2,0.2]]; %5 p = [p;[2.5,1.3]]; %6 p = [p;[2.2,2]]; %7 p = [p;[1.1,2]]; %8 % bloky TP = []; TP = [TP;[1 2 3 4]]; TP = [TP;[2 5 6 3]]; TP = [TP;[3 6 7 8]]; nx1 = 10; nx2 = 14; ny1 = 16; ny2 = 14; % deleni bloku N = []; N = [N;[nx1,ny1]]; N = [N;[nx2,ny1]]; N = [N;[nx2,ny2]]; % zahusteni bloku z = []; z = [z;[1,1,1,1]]; z = [z;[2,1,2,1]]; z = [z;[2,2,2,2]];
Obr. 1 Zobrazení jednotlivých bloků tvořících 2D řez sítě.

Pokud jsme s vlastní 2D sítí spokojeni, můžeme přejít k tvorbě 3D sítě. Tato síť se generuje příkazem generate3Dmesh.m, kde je nutné pouze stanovit počet buněk v obvodovém směru a excentricitu rotoru. Výsledná síť může vypadat např. takto, viz obrázek 2.

Obr. 2 2D a 3D výpočtová síť pro případ konkrétní ucpávky.

Vytvořená 3D výpočtová síť se uloží do adresáře mesh. Podobně 2D výpočtová síť se uloží do adresáře mesh2D. Formát výpočetních sítí v obou adresářích odpovídá formátu pro výpočtový software FlowPro. Vlastní výpočet probíhá nejprve na 2D síti. Po skončení výpočtu se z 2D výsledků extrapoluje 3D počáteční podmínka a provede se 3D výpočet. Pro extrapolaci dat z 2D do 3D se použije funkce to3D.m.

Inicializační 2D výpočet

Jak již bylo řečeno v předchozí sekci, k inicializaci počáteční podmínky ve 3D se využije 2D výpočet na 2D síti. K tomuto výpočtu lze použít standardní modul softwaru FlowPro, který se jmenuje Aerodynamics. Tento modul implementuje Navierovy-Stokesovy rovnice popisující laminární proudění stlačitelné vazké tekutiny a také několik turbulentních modelů. V případě použití Spalartova-Allmarasova modelu je třeba specifikovat následující parametry

model = 'SpalartAllmaras' dimension = %%% either kappa, or cp and cv must be defined %cp = % heat capacity at constant pressure cv = % heat capacity at constant volume kappa = % heat capacity ratio isFlowViscous = % true / false %%% for viscous flow (isFlowViscous = true) either reynolds and prandtl, %%% or viscosity and conductivity must be defined reynolds = % Reynolds number prandtl = 0.72 % Prandtl number %viscosity = % dynamic viscosity %conductivity = % thermal conductivity %%% inlet isInletSupersonic = % true / false %attackAngle = % angle of attack %%% subsonic inlet boundary condition (isInletSupersonic = false) pIn0 = % inlet stagnation pressure rhoIn0 = % inlet stagnation density %TIn0 = % inlet stagnation temperature (can be defined instead of rhoIn0) vtIn = % turbulence at the inlet %%% supersonic inlet boundary condition (isInletSupersonic = true) pIn = % inlet pressure rhoIn = % inlet density %TIn = % inlet temperature (can be defined instead of rhoIn) vIn = % inlet velocity %%% outlet pOut = % outlet pressure %machInf = % Mach number of undisturbed flow (can be defined instead of pOut) %%% turbulence model constant sigma = 0.6667 cb1 = 0.1355 cb2 = 0.622 ka = 0.41 cw2 = 0.3 cw3 = 2 cv1 = 7.1 ct1 = 1 ct2 = 2 ct3 = 1.2 ct4 = 0.5 Prt = 0.9 C_prod = 2

3D výpočet

Pro 3D výpočet je určený vytvořený modul ALabyrinthSeal3DA, který implementuje Navierovy-Stokesovy rovnice v rotujícím souřadnicovém systému. Pro modelování turbulentní viskozity je zde uvažován Spalartův-Allmarasův turbulentní model. Před vlastním výpočtem je třeba nadefinovat následující parametry modelu.

model = 'SpalartAllmarasNavierStokes3DRotFrame' dimension = 3 %%% either kappa, or cp and cv must be defined cp = % heat capacity at constant pressure cv = % heat capacity at constant volume % kappa = % heat capacity ratio %tEnd = % physical time at which the computatin is to be terminated isFlowViscous = true % true / false %%% for viscous flow (isFlowViscous = true) either reynolds and prandtl, %%% or viscosity and conductivity must be defined %reynolds = % Reynolds number %prandtl = % Prandtl number viscosity = % dynamic viscosity conductivity = % thermal conductivity %%% inlet isInletSupersonic = false % true / false attackAngleAlfa = % angle of attack attackAngleBeta = % angle of attack %%% subsonic inlet boundary condition (isInletSupersonic = false) pIn0 = % inlet stagnation pressure %rhoIn0 = % inlet stagnation density TIn0 = % inlet stagnation temperature (can be defined instead of rhoIn0) vtIn = %%% supersonic inlet boundary condition (isInletSupersonic = true) pIn = % inlet pressure rhoIn = % inlet density %TIn = % inlet temperature (can be defined instead of rhoIn) vIn = % inlet velocity lRef = %%% outlet pOut = % outlet pressure %machInf = % Mach number of undisturbed flow (can be defined instead of pOut) %%% turbulence model constant sigma = 0.6667 cb1 = 0.1355 cb2 = 0.622 ka = 0.41 cw2 = 0.3 cw3 = 2 cv1 = 7.1 ct1 = 1 ct3 = 1.2 ct2 = 2 ct4 = 0.5 Prt = 0.9 C_prod = 2 %%% omega omega = 0, 0, 0 omegaRotor = 0, 0, 0 omegaInlet = 0, 0, 0 rotorExcentricity = 0, 0, 0

Vyhodnocení dat

Software FlowPro standardně provádí export výsledků do txt souborů, popř. do vtk souboru. Při použití softwaru na vizualizaci vtk souborů (např. paraview) lze zobrazit rychlostní a tlakové pole. Klíčovou roli v labyrintových ucpávkách hraje obvodová změna tlakového pole v jednotlivých kavitách a také celkové síly působící na rotor. Z tohoto důvodu byly v systému Matlab implementovány následující skripty pro zpracování dat:

Příklad

Jako příklad si ukážeme výpočet ucpávky, jejíž geometrie je zobrazena na obr. \ref{mesh}. Jedná se vicekomorovou (labyrintovou) VT rotorovou ucpávku. Bloková výpočtová síť obsahuje 2404 buněk v řezu a 90 po obvodu, což odpovídá 216 360 buňkám a 865 440 stupňům volnosti (DGFEM druhého řádu). Pro výpočet 2D úlohy byly nastaveny parametry odpovídající proudění vzduchu s tlakovým spádem 6bar.

model = 'NavierStokes' dimension = 2 %%% either kappa, or cp and cv must be defined cp = 1005 % heat capacity at constant pressure cv = 717.9 % heat capacity at constant volume % kappa = 1.4 % heat capacity ratio isFlowViscous = true % true / false %%% for viscous flow (isFlowViscous = true) either reynolds and prandtl, %%% or viscosity and conductivity must be defined %reynolds = % Reynolds number %prandtl = 0.72 % Prandtl number viscosity = 1.879e-5 % dynamic viscosity conductivity = 0.0242 % thermal conductivity %%% inlet isInletSupersonic = false % true / false attackAngle = 0 % angle of attack %%% subsonic inlet boundary condition (isInletSupersonic = false) pIn0 = 13e5 % inlet stagnation pressure %rhoIn0 = % inlet stagnation density TIn0 = 295.3 % inlet stagnation temperature (can be defined instead of rhoIn0) %%% supersonic inlet boundary condition (isInletSupersonic = true) pIn = % inlet pressure rhoIn = % inlet density %TIn = % inlet temperature (can be defined instead of rhoIn) vIn = % inlet velocity lRef = 1e-3 %%% outlet pOut = 715000 % outlet pressure %machInf = % Mach number of undisturbed flow (can be defined instead of pOut)
Obr. 3 Machovo číslo a tlakové pole pro inicializační 2D výpočet.

Na obr. 3 je znázorněno rozložení Machova čísla a tlakového pole v ucpávce spočtené pomocí 2D simulace. Použitím funkce to3D.m převedeme 2D data do 3D. Výsledkem je soubor initW.txt, který představuje počáteční podmínku ve 3D. 3D výpočet lze provést přímo, bez generování počáteční podmínky, nicméně takto generovaná 3D počáteční podmínka poměrně výrazně zkracuje dobu potřebnou pro konvergenci výpočtu. Pro 3D výpočet jsme volili následující nastavení parametrů, kde jsme uvažovali frekvenci precese 20Hz, otáčky rotoru 7496 ot/min a excentricitu e=3e-5 m.

model = 'SpalartAllmarasNavierStokes3DRotFrame' dimension = 3 %%% either kappa, or cp and cv must be defined cp = 1005 % heat capacity at constant pressure cv = 717.9 % heat capacity at constant volume % kappa = 1.4 % heat capacity ratio %tEnd = % physical time at which the computatin is to be terminated isFlowViscous = true % true / false %%% for viscous flow (isFlowViscous = true) either reynolds and prandtl, %%% or viscosity and conductivity must be defined %reynolds = % Reynolds number %prandtl = 0.72 % Prandtl number viscosity = 1.879e-5 % dynamic viscosity conductivity = 0.0242 % thermal conductivity %%% inlet isInletSupersonic = false % true / false attackAngleAlfa = 0 % angle of attack attackAngleBeta = 0 % angle of attack %%% subsonic inlet boundary condition (isInletSupersonic = false) pIn0 = 13e5 % inlet stagnation pressure %rhoIn0 = % inlet stagnation density TIn0 = 295.3 % inlet stagnation temperature (can be defined instead of rhoIn0) vtIn = 0.01 %%% supersonic inlet boundary condition (isInletSupersonic = true) pIn = % inlet pressure rhoIn = % inlet density %TIn = % inlet temperature (can be defined instead of rhoIn) vIn = % inlet velocity lRef = 1e-3 %%% outlet pOut = 715000 % outlet pressure %machInf = % Mach number of undisturbed flow (can be defined instead of pOut) %%% turbulence model constant sigma = 0.6667 cb1 = 0.1355 cb2 = 0.622 ka = 0.41 cw2 = 0.3 cw3 = 2 cv1 = 7.1 ct1 = 1 ct3 = 1.2 ct2 = 2 ct4 = 0.5 Prt = 0.9 C_prod = 2 %%% omega omega = 125.6, 0, 0 omegaRotor = 785, 0, 0 omegaInlet = 392.5, 0, 0 rotorExcentricity = 0, 0, 0.03

Na obr. 4 je znázorněno rozložení Machova čísla a tlakového pole v ucpávce spočtené pomocí 3D simulace a vizualizované pomocí softwaru Paraview. Nyní můžeme pomocí skriptů v Matlabu analyzovat výsledky. Amplitudy tlaku a rychlosti jsou znázorněny na obrázcích 5 a 6. Swirl rychlosti je zobrazený na obr. 7. Na obr. 8 je znázorněna amplituda tlaku na rozvinuté válcové ploše, jejíž osa symetrie je totožná s osou statoru. V tabulce \ref{tab} jsou uvedeny hodnoty amplitudy tlaku v předem definovaných bodech o souřadnicích [xj,yj]. Nejdůležitější funkcí je funkce computeForce.m, která z tlakového pole vypočítá sílu jakou proudící tekutina působí na rotor. Voláním této funkce dostaneme F [N] = [1056.3, 4.0, -12.7]. Provedením několika simulací pro různé hodnoty precese je možné ze získaných silových účinků vypočítat hodnoty koeficientů tuhosti a tlumení ucpávky, které je dále možné využít při simulacích dynamiky rotoru.

Obr. 4 Machovo číslo a tlakové pole spočtené 3D simulací.
Obr. 5 pressure2DamplitudeOnSurface.m - amplituda tlaku zobrazena na válcové ploše s osou rotace totožnou s osou statoru.
Obr. 6 pressure3Dto2Damplitude.m - amplituda tlaku v logaritmické stupnici.
Obr. 7 velocity3Dto2Damplitude.m - amplituda rychlosti v logaritmické stupnici.
Obr. 8 velocitySwirl.m - swirl rychlosti.