Mathematical Physics - Volume II - Numerical Methods

Chapter 4. Finite volume method

118

the system (4.13)-(4.15), supplemented by the relation (4.5) is set up. With known values of the dependent variable q n inside each cell ( i , j , k ) , approximate form of the system of equations (4.13.1) is obtained by I + β ∆ t ( δ ξ A n + δ η B n + δ ζ C n ) ∆ q n + ∆ t R n = 0 , (4.19) where the residual R n is R n = δ ξ F ( q n )+ δ η G ( q n )+ δ ζ H ( q n ) . (4.19.1) In the equations (4.19) and (4.19.1) δ ξ , δ η and δ ζ represent central difference operators ∂ / ∂ ξ , ∂ / ∂η and ∂ / ∂ ζ . The parameter β in the expression (4.19) determines time accuracy of applied scheme. For the value β = 0 . 5 the scheme is second order accurate, while accuracy decreases to the first order for β = 1. Unfactored implicit scheme, based on the equation (4.19), requires large available computer memory space for solving huge band block matrices, which is very difficult to provide for the case of three-dimensional flows. Unconditionally stable implicit scheme, with the largest possible error up to the order ( ∆ t ) 2 , regardless of the number of spatial dimensions, is derived by the so-called LU factorization [20], [21] and [6] I + β ∆ t ( δ − ξ A + + δ − η B + + δ − ζ C + ) n ∗ ∗ I + β ∆ t ( δ + ξ A − + δ + η B − + δ + ζ C − ) n ∆ q n + ∆ t R n = 0 , (4.20) where δ − ξ , δ − η and δ − ζ are difference operators “backward” oriented, while δ + ξ , δ + η and δ + ζ are “forward” oriented: δ − ξ ( A + ∆ q n ) i , j , k = A + i + 1 / 2 , j , k ∆ q n i , j , k − A + i − 1 / 2 , j , k ∆ q n i − 1 , j , k , δ + ξ ( A − ∆ q n ) i , j , k = A − i + 1 / 2 , j , k ∆ q n i + 1 , j , k − A − i − 1 / 2 , j , k ∆ q n i , j , k , (4.20.1) where the terms of the matrices with indices ( i + 1 / 2 , j , k ) and ( i − 1 / 2 , j , k ) are computed by averaging flow variables between grid cells defined by indices ( i , j , k ) and ( i + 1 , j , k ) , ie. ( i − 1 , j , k ) and ( i , j , k ) . Difference operators related to the remaining two spatial directions can be determined on the basis of similar relations. Matrix decomposition in factorized scheme (4.20) provides diagonal dominance of terms within the parentheses, while implicitly introduces an artificial viscosity, necessary to stabilize the applied scheme based on central difference operators. Flux matrices A + , B + , C + , A − , B − and C − in the equation (4.20) are formed in way providing that eigenvalues of the “+” matrix are non-negative, and the eigenvalues of “-” matrix are non-positive

1 2

1 2

A − =

A +

( A + r A I ) ,

( A − r A I ) , ( B − r B I ) , ( C − r C I ) ,

=

1 2

1 2

B − =

B +

(4.20.2)

( B + r B I ) ,

=

1 2

1 2

C − =

C +

( C + r C I ) ,

=

where I is the unit matrix, while the factors r A , r B and r C are defined as follows: r A ≥ max ( | λ A | ) , r B ≥ max ( | λ B | ) and r C ≥ max ( | λ C | ) ,

(4.20.3)

bearing in mind that λ A , λ B and λ C represent eigenvalues of flux matrices A , B and C , respectively,