方程及状态参数
进排气管内一维变截面非定常控制方程为:
{ ∂ ( ρ A ) ∂ t + ∂ ( ρ u A ) ∂ x = 0 ∂ ( ρ u A ) ∂ t + ∂ ( ρ u 2 A + p A ) ∂ x = − ρ G A + p ∂ A ∂ x ∂ ( ρ e 0 A ) ∂ t + ∂ ( ρ u h 0 A ) ∂ x = ρ q ˙ A (1) \left\{ \begin{array}{l} \frac{
{\partial \left( {\rho A} \right)}}{
{\partial t}} + \frac{
{\partial \left( {\rho uA} \right)}}{
{\partial x}} = 0\\ \frac{
{\partial \left( {\rho uA} \right)}}{
{\partial t}} + \frac{
{\partial \left( {\rho {u^2}A + pA} \right)}}{
{\partial x}} = - \rho GA + p\frac{
{\partial A}}{
{\partial x}}\\ \frac{
{\partial \left( {\rho {e_0}A} \right)}}{
{\partial t}} + \frac{
{\partial \left( {\rho u{h_0}A} \right)}}{
{\partial x}} = \rho \dot qA \end{array} \right. \tag{1} ⎩⎪⎨⎪⎧∂t∂(ρA)+∂x∂(ρuA)=0∂t∂(ρuA)+∂x∂(ρu2A+pA)=−ρGA+p∂x∂A∂t∂(ρe0A)+∂x∂(ρuh0A)=ρq˙A(1)
其中总热力学能为静热力学能和动能之和: e 0 = e + 1 2 u 2 = R g T γ − 1 + 1 2 u 2 = p ρ ( γ − 1 ) + 1 2 u 2 {e_0} = e + \frac{1}{2}{u^2} = \frac{
{
{R_g}T}}{
{\gamma - 1}} + \frac{1}{2}{u^2} = \frac{p}{
{\rho \left( {\gamma - 1} \right)}} + \frac{1}{2}{u^2} e0=e+21u2=γ−1RgT+21u2=ρ(γ−1)p+21u2.
令 U = [ ρ A ρ u A ρ e 0 A ] {\bf{U}} = \left[ \begin{array}{l} \rho A\\ \rho uA\\ \rho {e_0}A \end{array} \right] U=⎣⎡ρAρuAρe0A⎦⎤, F = [ ρ u A ( ρ u 2 + p ) A ρ h 0 u A ] {\bf{F}} = \left[ \begin{array}{l} \rho uA\\ \left( {\rho {u^2} + p} \right)A\\ \rho {h_0}uA \end{array} \right] F=⎣⎡ρuA(ρu2+p)Aρh0uA⎦⎤, S = [ 0 p d A d x − ρ G A ρ q ˙ A ] {\bf{S}} = \left[ \begin{array}{l} 0\\ p\frac{
{dA}}{
{dx}} - \rho GA\\ \rho \dot qA \end{array} \right] S=⎣⎡0pdxdA−ρGAρq˙A⎦⎤,则式(1)化为:
∂ U ∂ t + ∂ F ( U ) ∂ x = S ( U ) (2) \frac{
{\partial {\bf{U}}}}{
{\partial t}} + \frac{
{\partial {\bf{F}}\left( {\bf{U}} \right)}}{
{\partial x}} = {\bf{S}}\left( {\bf{U}} \right)\tag{2} ∂t∂U+∂x∂F(U)=S(U)(2)
在初始化网格热力状态时,一般给的都是 u , p , T u,p,T u,p,T和广义组分表征量 α \alpha α,然后求 U \bf{U} U,这里定义一个类Node,Node中的一个方法为 U {\bf{U}} U的初始化函数:
def Uinit(self, u, p, T, AFAK=1.e8, A=1):
from GasProperty import Rg
from GasProperty import k_Justi
self.p = p
self.T = T
self.rho = p / (Rg(AFAK) * T)
self.u = u
self.k = k_Justi(self.T, AFAK)
self.e = self.p / self.rho / (self.k - 1.)
self.e0 = self.e + self.u ** 2 / 2.
self.h0 = self.e0 + self.p / self.rho
self.a = pow(self.k * self.p / self.rho, 1. / 2.)
from numpy import array
self.U = array([self.rho * A, self.rho * self.u * A,
(self.p / (self.k - 1) + 1. / 2. * self.rho * self.u ** 2) * A]).transpose()
F ( U ) {\bf{F}}(\bf{U}) F(U)的表达式:
F ( U ) = [ U 2 U 2 2 U 1 + ( γ − 1 ) ( U 3 − U 2 2 2 U 1 ) U 2 ( U 3 + ( γ − 1 ) ( U 3 − U 2 2 2 U 1 ) ) U 1 ] (3) {\bf{F}({\bf{U}})} = \left[ \begin{array}{l} {U_2}\\ \frac{
{
{U_2}^2}}{
{
{U_1}}} + \left( {\gamma - 1} \right)\left( {
{U_3} - \frac{
{
{U_2}^2}}{
{2{U_1}}}} \right)\\ \frac{
{
{U_2}\left( {
{U_3} + \left( {\gamma - 1} \right)\left( {
{U_3} - \frac{
{
{U_2}^2}}{
{2{U_1}}}} \right)} \right)}}{
{
{U_1}}} \end{array} \right]\tag{3} F(U)=⎣⎢⎢⎢⎡U2U1U22+(γ−1)(U3−2U1U22</