Mathematical Physics - Volume II - Numerical Methods

4.2 Solution of Euler equations

113

Special attention has to be paid to the analysis of dissipative terms behavior at the boundary of the physical domain. Taking into account relations (4.9.5) and (4.9.6), introduction of second order terms requires information of the value of flow variables in neighboring grid cells, while in the case of fourth-order terms influence domain extends to two adjacent cells in each coordinate direction. For cells at the boundary of the physical domain defining of the artificial viscosity terms is achieved by combining physical and numerical boundary conditions. In this way the lack of necessary information is compensated, and at the same time the stability of the scheme at the boundaries is ensured. After the introduction of dissipative terms the system of equations (4.7) receives the final form ∂ ∂ t ( h i , j , k W i , j , k )+ Q ( W ) i , j , k − D ( W ) i , j , k = 0 . (4.10) Since the volume of the hexahedral cell is invariant in time, application of Runge-Kutta method for solving system of ordinary differential equations leads to relations W ( 0 ) i , j , k = W n i , j , k , W ( r ) i , j , k = W ( 0 ) i , j , k − α r h − 1 i , j , k ∆ t [ Q ( W ) ( r − 1 ) i , j , k − D ( W ) ( r − 1 ) i , j , k ] , r = 1 , 2 , . . . , r max W n + 1 i , j , k = W ( r max ) i , j , k , (4.11) where r max determines the order of the adopted Runge-Kutta scheme, while α r are the corresponding coefficients of the scheme. In equations (4.11) the indices n and n + 1 indicate the values of the variable W i , j , k at the beginning and the end of an iterative cycle. To achieve significant savings in computer time, artificial viscosity terms are most often "frozen" on the first, ie. the second step within an iterative cycle, bearing in mind that the stability of such a modified equation is almost unchanged. If the calculation of dissipative terms is performed only within the first step of each iterative cycle, the system of equations (4.11) receives form W ( 0 ) i , j , k = W n i , j , k , W ( r ) i , j , k = W ( 0 ) i , j , k − α r h − 1 i , j , k ∆ t [ Q ( W ) ( r − 1 ) i , j , k − D ( W ) ( 0 ) i , j , k ] , r = 1 , 2 , . . . , r max W n + 1 i , j , k = W ( r max ) i , j , k . (4.11.1) When the values inside the brackets on the right hand of equation (4.11) and equation (4.11.1), become zero, a steady state is reached and the iterative process is completed. The steady state does not depend on the value of the adopted integration step ∆ t . In this case, the value of the variable W ( r max ) i , j , k at the end of iterative cycle is equal to the value of W ( 0 ) i , j , k at the beginning of the iteration, ie W n + 1 i , j , k = W n i , j , k . A very important detail in this procedure is determination of the integration step ∆ t , taking into account that it is the largest value determined by the time interval required to perturbation transfers from one cell boundary to another. Integration step ∆ t for the cell ( i , j , k ) is determined by the relation ∆ t i , j , k = 1 ( ∆ t ξ ) i , j , k + 1 ( ∆ t η ) i , j , k + 1 ( ∆ t ζ ) i , j , k − 1 , (4.12) in which ( ∆ t ξ ) i , j , k , ( ∆ t η ) i , j , k and ( ∆ t ζ ) i , j , k are time intervals of disturbance propagation in the direc tions of the corresponding coordinate axes within the boundary of the observed grid cell. By mapping to a curvilinear coordinate system ( ξ , η , ζ ) , in which constant values of curvilinear coordinates