Aeroblade is a set of subroutines designated for studying of fluid flow through a blade cascade in 2D. The subroutines consist of blade dynamic models implemented in FlowPro and MATLAB scripts, which is used for data processing. For download click here.

Úvod

Softwarový modul "Aeroblade 2D" 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 "Aeroblade 2D" umožňuje analýzu aeroelastických jevů a posouzení stability olopatkovaného disku při samobuzeném kmitání ve 2D řezu. Výpočtový modul obsahuje rovněž programové nástroje pro vyhodnocení vypočtených dat, např. pro stanovení tlumícího koeficientu.

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

Před vlastním výpočtem je vždy nutné vygenerovat výpočetní síť. Tato činnost je však vždy velmi časově náročná. Z tohoto důvodu byl v rámci modulu "Aeroblade 2D" vyvinut automatický generátor sítě. Generátor je naprogramovaný jako skript v jazyce Matlab. Do generátoru vstupují tři základní parametry sítě, tj. body profilu, posunutí profilu H a počet profilů n. Na obrázku 1 je znázorněna výpočtová síť okolo 3 profilů SE1050 s posunutím H = 0.3. Kromě výpočtové sítě se dále automaticky generují blending funkce (obrázek 2), které jsou třeba pro výpočet deformované sítě.

Obr. 1 Automaticky generovaná čtyřúhelníková síť s profilu SE10510, s posunutím H = 0.3 a spočtem 3 profilů. Černá tečka vyznačuje střed hmotnosti profilů (vlevo). Body profilu SE1050 (vpravo).

Obr. 2 Automaticky generované blending funkce okolo jednotlivých profilů, soužící pro přepočet sítě.

Nastavení vstupních parametrů výpočtu

Pomocí FlowPro manageru vybereme modul Aerodynamics, který implementuje základní rovnice proudění. Parametry výpočtového modelu jsou specifikované v souboru parameters.txt. Pro systém Navierových-Stokesových rovnic vypadají parametry modelu následovně:

model = 'NavierStokes' dimension = 2 %%% 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 = % 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) %%% 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)

V nabídce FlowPro manageru dále zvolíme modul Rigid2D, který implementuje lineární dynamický model tuhého tělesa (lopatky). Parametry tohoto modelu vypadají následovně:

dynamicsModel = 'Rigid2D' meshDeformationType = 'Rigid2D' %zLength = % length in z direction %tKick = % dimensionalless time of kinematic motion (default is inf) %rotationCenters = % coordinates of gravity centers %M = % mass matrixes %B = % damping matrixes %K = % stiffness matrixes %xMotion = % kinematic motion in x direction %yMotion = % kinematic motion in y direction %alphaMotion = % kinematic motion in alpha direction %xExternForceInTime = % time force in x direction %yExternForceInTime = % time force in y direction %alphaExternForceInTime = % time force in alpha direction

Testovací úloha

V této sekci je stručně uvedeno několik testovacích úloh modulu "Aeroblade 2D", které validují vytvořený modul pro úlohy interakce lopatek s proudící tekutinou.

Kinematicky buzený profil NACA0012

V této testovací úloze sledujeme účinky proudění stlačitelné nevazké tekutiny na kmitající profil NACA0012. Kmitající profil je umístěn do nerozrušeného proudu definovaného Machovým číslem M_∞ = 0.755. Časový průběh úhlu natočení profilu je popsán funkcí

\alpha(t) = 0.016^{\circ} + 2.51^{\circ}\sin(\omega\,t),

kde úhlová rychlost je definována jako ω = 2ku_∞/c, kde k = 0.0814 je redukovaná frekvence, délka tětivy je c = 1 a u_∞ je rychlost nerozrušeného proudu. Na obrázku 3 (vpravo) je vidět dobrá shoda mezi vypočteným vztlakovým koeficientem a hodnotami z literatury [D.J. Kirshman, F. Liu: Computers and Fluids, 35(6): 571-586, 2006].

Obr. 3 Profil NACA0012 (vlevo), vztlakový koeficient v závislosti na úhlu natočení (vpravo).

Hranice flutteru

V této testovací úloze je rekonstruována hranice flutteru pro letecký profil NACA 64A-010. Dynamický model lopatky, viz obrázek 4, byl převzat z literatury [K. Isogai:AIAA Journal, 17 (7), pp. 793-795].

\mathbf{M}\ddot{\bsym{y}} + \mathbf{K}\bsym{y} = \mathbf{F}.

Prvky matice hmotnosti a tuhosti jsou definovány jako

\mathbf{M} = \left[ \begin{array}{cc} m & S_{\alpha} \\ S_{\alpha} & I_{\alpha} \end{array} \right], \hspace{6pt} \mathbf{K} = \left[ \begin{array}{cc} K_h & 0 \\ 0 & K_{\alpha} \end{array} \right],\hspace{3pt}\mathbf{B} = \mathbf{0},\hspace{3pt} \mathbf{F} = [L,M_A]^T,\hspace{3pt} \mathbf{y} = [h,\alpha]^T,

m = \mu \pi \rho_{\infty} b^2,\hspace{2pt} I_{\alpha} = m b^2 r_{\alpha}^2,\hspace{2pt} S_{\alpha} = m b x_{\alpha}, \hspace{2pt} \omega_{\alpha} = \frac{U_{\infty}}{U_F b \sqrt{\mu}},\hspace{2pt} K_h = m\omega_h^2, \hspace{2pt} K_{\alpha} = I_{\alpha}\omega_{\alpha}^2,

s parametry

a = -2,\hspace{12pt} b = 0.5,\hspace{12pt} x_{\alpha} = 1.8,\hspace{12pt} r_{\alpha}^2 = 3.48,\hspace{12pt} \omega_{\alpha} = \omega_{h},\hspace{12pt} \mu = 60, \nonumber

kde U_F je rychlostní index a ρ_∞, U_∞ jsou hodnoty hustoty a rychlosti nerozrušeného proudu, které jsou funkcí Machova čísla.

Obr. 4 Isogaiův dynamický model lopatky.

Hranice flutteru je definována jako relace mezi rychlostním indexem U_F a Machovým číslem M_∞ taková, že pro každou dvojici [U_F, M_∞] nedochází ani k nestabilitě typu flutter ani k tlumení počátečních výchylek profilu. Tato hranice se určuje iteračním postupem. Obrázek 5 ukazuje časový vývoj výchylek pro Machovo číslo M_∞=0.75 a rychlostní index U_F = 1 (vlevo), kde dochází k tlumení výchylek a U_F = 1.25 (vpravo), kde dochází k nestabilitě typu flutter. Hranice flutteru bude ležet někde mezi těmito rychlostními indexy. Na obrázku 6 (vlevo) je rekonstruovaná hranice flutteru pro NACA 64A-010 s Isogaiovým dynamickým modelem lopatky. Porovnání s literaturou ukazuje velmi dobrou shodu.

Obr. 5 Časový vývoj výchylky $h$ a úhlu natočení $\alpha$ v závislosti na Machově čísle M_∞=0.75 a rychlostním indexu U_F = 1 (vlevo) a U_F = 1.25 (vpravo).

Obr. 6 Rekonstrukce hranice flutteru a porovnání s literaturou (vlevo), síť okolo NACA 64A-010 profilu (vpravo).

Výpočet, zpracování a vizualizace výsledků

Vlastní výpočet mapy flutteru a hledání hranice flutteru je poměrně složitý výpočet. Je třeba provádět velké množství simulací, které je třeba průběžně vyhodnocovat a podle průběžných výsledků dále upravovat parametry dalších simulací. Z tohoto důvodu byl v Matlabu vytvořený skript, který tyto simulace provádí automaticky. Vstupem do skriptu je vektor Machových čísel M_∞,k, k = 1,...,n a hraniční hodnoty indexu rychlosti [U_F,min,U_F,max]. Algoritmus pak pro každou hodnotu Machova čísla M_∞,k pomocí metody půlení intervalu stanový hodnotu rychlostního indexu U_F,k na hranici flutteru.

Obr. 7 Testovací profil ploché lopatky (vlevo), automaticky generovaná síť okolo testovacích profilů (vpravo).

Jako příklad uvažujme profil na obrázku 7 (vlevo). Pro generování síťe profil natočíme o 20°, počet profilů zvolíme 4 a posunutí volíme H=0.3. Výpočetní síť generovaná automatickým generátorem je vidět na obrázku 7 (vpravo). Parametry výpočtu se nastavují automaticky před každým výpočtem. Soubor parameters.txt může vypadat následovně.

%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% physical parameters %%% %%%%%%%%%%%%%%%%%%%%%%%%%%% model = 'NavierStokes' dimension = 2 %%% exactly two out of the three following parameters must be defined: kappa, cp, cv %cp = % heat capacity at constant pressure cv = 1 % heat capacity at constant volume kappa = 1.4 % heat capacity ratio isFlowViscous = false % 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 = false % true / false attackAngle = -0.3491 % angle of attack %%% subsonic inlet boundary condition (isInletSupersonic = false) pIn0 = 1 % inlet stagnation pressure rhoIn0 = 1 % inlet stagnation density %TIn0 = % 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 %%% outlet %pOut = 0.9 % outlet pressure machInf = 0.95 %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% numerical parameters %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% continueComputation = false % true / false % terminating conditions steps = 3000 % computation is terminated after defined number of time steps is reached %residuum = % computation is terminated after defined residuum is reached (optional) %tEnd = % computation is terminated after defined physical time is reached (optional) saveRate = 10 % save rate animation = true % true/false movingMesh = true % true for ALE computation order = 2 % order of accuracy in space orderInTime = 2 % order of accuracy in time (1 or 2 is supported) CFL = 5 threads = 6 % number of CPU threads that will be used for computation newtonIters = 1 penalty = 1 % damping dampTol = 3 dampConst = 0 % parallel mode overlap = 2 % must be specified only for parallel computations schwarzIters = 0 schwarzTol = 1e-10 %%%%%%%%%%%%%%%%%%%%%%%%%%% %%% bodies parameters %%% %%%%%%%%%%%%%%%%%%%%%%%%%%% dynamicsModel = 'Rigid2DGACR' meshDeformationType = 'Rigid2D' simID = M0p95UF0p175 %zLength = % length in z direction %tKick = 1 % dimensionalless time of kinematic motion (default is inf) tKickForce = 3 rotationCenters = 0, 0.15374, 0, 0.45374, 0, 0.75374, 0, 1.0537 % gravity centres M = 1, 0, 0, 0, 10.3739, 5.187, 0, 5.187, 9.0253 B = 0, 0, 0, 0, 0, 0, 0, 0, 0 K = 0, 0, 0, 0, 72.5112, 0, 0, 0, 63.0848 Kint = 0, 0, 0, 0, 36.2556, 0, 0, 0, 0 %xMotion = % kinematic motion in x direction %yMotion = % kinematic motion in y direction %alphaMotion = % kinematic motion in alpha direction %xExternForceInTime = % time force in x direction yExternForceInTime = 0.001*Math.sin(2.6438*t),0.001*Math.sin(2.6438*t + Math.PI/2), 0.001*Math.sin(2.6438*t + 2*Math.PI/2),0.001*Math.sin(2.6438*t + 3*Math.PI/2) alphaExternForceInTime = 0.001*Math.sin(2.6438*t),0.001*Math.sin(2.6438*t + Math.PI/2), 0.001*Math.sin(2.6438*t + 2*Math.PI/2),0.001*Math.sin(2.6438*t + 3*Math.PI/2)

Vlasní výpočet se spouští v Matlabu příkazem computeBatchSplitInterval. Po ukončení výpočtu lze zobrazit výsledky příkazem processData. Pro uvažovanou úlohu jsou výsledky zobrazeny na obrázku 8.

Obr. 8 Hranice flutteru pro testovanou kaskádu. Červené tečky značí nestabilní oblast a zelené stabilní.