Mathematical Physics - Volume II - Numerical Methods

4.3 Solution of Navier–Stokes equations

151

where Jacobi transformation matrix J = ∂ ( x , y , z ) / ∂ ( ξ , η , ζ ) is calculated from the expres sion (4.15.6). For the perfect gas system of equations (4.126) is completed by defining fluid energy (4.5), thus finally closing the Navier-Stokes system of equations. In the system of equations (4.126), ie. (4.126.1) the unknown to be determined is a vector of variable flow quantities q . After the physical flow domain is discretized, as described in the chapter 4.2.2, a system of equations that approximates partial differential system of equation (4.126), can be set up. If inside each computational cell ( i , j , k ) the value of the dependent variable q n is known, approximate form of the system of equations (4.126.1) can be obtained n I + β ∆ t h δ ξ A n + δ η B n + δ ζ ( C n − C n v ) io ∆ q n + ∆ t R n = 0 , (4.128) where the matrices A , B and C are previously defined by expressions (4.2.2), while the residual term is R n R n = δ ξ F ( q n )+ δ η G ( q n )+ δ ζ H ( q n ) − H v ( q n ) . (4.128.1) In the equations (4.128) and (4.128.1) δ ξ , δ η and δ ζ represent central difference opera tors ∂ / ∂ ξ , ∂ / ∂η and ∂ / ∂ ζ . The parameter β in the expression (4.128) determines time accuracy of the applied scheme, which is analyzed in chapter 4.2.2. Unconditionally stable implicit scheme is derived by LU factorization [20], [21] and [6] n I + β ∆ t h δ − ξ A + + δ − η B + + δ − ζ ( C + − C + v ) io n ∗ ∗ n I + β ∆ t h δ + ξ A − + δ + η B − + δ + ζ ( C − − C − v ) io n ∆ q n + ∆ t R n = 0 , (4.129) where δ − ξ , δ − η , δ − ζ , δ + ξ , δ + η and δ + ζ are difference operators determined by relations (4.20.1), where the elements of the matrices with indices ( i + 1 / 2 , j , k ) and ( i − 1 / 2 , j , k ) are cal culated by averaging flow variables between computational cells ( i , j , k ) and ( i + 1 , j , k ) , ie. ( i − 1 , j , k ) and ( i , j , k ) , respectively. Difference operators related to the remaining two coordinate directions can be determined on the basis of similar relations. Solutions of the system of equations (4.129) is possible to determine in two passes: n I + β ∆ t h δ − ξ A + + δ − η B + + δ − ζ ( C + − C + v ) io n ∆ q ∗ n = − ∆ t R n , n I + β ∆ t h δ + ξ A − + δ + η B − + δ + ζ ( C − − C − v ) io n ∆ q n = ∆ q ∗ n , (4.129.1) where the matrices of “viscous” flux C + v and C − v , present in the implicit term of the equation (4.129), are approximated according to the work of Pulliam [37] C + v = λ C v I i C − v = − λ C v I , (4.129.2) where λ C v is the eigenvalue of the viscous flux matrix C v λ C v = µ ( ζ 2 x + ζ 2 y + ζ 2 z ) ∂ ζ 1 ρ , (4.129.3)