Implicit boundary equations for conservative Navier–Stokes equations
Introduction
Modern CFD simulations solve the Navier–Stokes (NS) equations on discretized computational grid points or cells within finite flow domains in a general curvilinear coordinate system: The definition of each term and some typical numerical procedures of time-dependent solutions, either explicit or implicit, can be found in textbooks [1]. In general, the stability constraint of explicit time integration schemes imposes a stringent limit on the maximum allowable time steps, and consequently results in excessive numbers of iterations and long CPU time for stiff problems or very fine computational grids. In comparison, implicit time integration schemes are more attractive, because they are usually unconditionally stable to allow large time steps, and the only constraint is the desired temporal resolution of the physical problems.
The implicit solution of Eq. (1) is challenging, and the solution techniques will be described in section 2. However, many issues can be illustrated using the scalar nonlinear Burger's equation: where, R is again referred to as the residuals involving spatial derivatives but in the x direction only, and the flux is . The numerical solutions are sought at discretized time instances and grid points within , subjecting to certain initial and boundary conditions.
Applying the first-order backward Euler method to this partial differential equation (PDE) leads to the following semi-discretized equation:
Eq. (3) can not be solved directly due to the nonlinear terms in the residuals R and the flux f. One way is to linearize these terms at the previous time instance using the Taylor series expansion, and organize the resulting equation into a delta form: where, and the numerical solutions will be . Or alternatively, one may re-write Eq. (3) into , and then solve it using the Newton iteration: where, the iteration can start with , and accurate solutions are recovered as when converges to machine zero.
Generally, numerical solutions of Eq. (4) suffer from linearization errors, and those of Eq. (5a), (5b) are expected to be accurate. Either way, one needs to solve the equation: where, the left implicit operator (in particular, the Jacobian matrix ) depends on the spatial discretization schemes. For example, when finite difference schemes are used, where, l and r indicate the left and right shifts of the grid points in the stencil. Therefore, and with , the left matrix can be readily assembled.
However, the equations above are not applicable on boundaries, because the stencils of the spatial discretization schemes protrude through the boundaries, and more importantly, a well-posed physical problem demands special boundary conditions [2]. For example, if in the entire domain , at the inlet must be given to establish the problem, and this is referred to as the physical boundary conditions. At the outlet , although no boundary conditions are need to solve the original PDE, a numerical procedure is necessary to determine the unknown boundary values, and this is referred to as the numerical boundary conditions. One may simply extrapolate the boundary value from interior solutions using, for instance, , or switch to one-sided scheme, which essentially regards Eq. (2) as the exact governing equation at the outlet. As a result, the governing equations and their discretized forms on boundaries are usually dramatically different from those above for interior grids.
Discussions above can be applied to the NS equations (1) directly, except two extra complications arise: (i) the increased dimensions of the solution vector and their spatial derivatives in other two directions inevitably yield more complicated left implicit operators. Direct inversion of the large difference equation systems appears difficult or even impossible. Many solution techniques have been developed to solve the discretized implicit equation systems, including the LU-SGS method [3], [4], the data parallel lower–upper relaxation (DPLUR) method [5] and the data parallel linear relaxation method [6], the Newton iteration combined with the GMRES method [7], [8], [9]; (ii) the original NS equations do not satisfy all boundary conditions automatically. Although the physical boundary conditions can be expressed in terms of some known boundary values, numerical boundary conditions cannot be set arbitrarily. Therefore, coexistence of the physical and numerical boundary conditions will make the boundary treatment a more challenging task. Ideally, boundary conditions in terms of the conservative variables are preferred, so that they can be discretized into and then collectively solved with those of the inner grids [10].
A survey of existing publications suggest that there are two major approaches of implicit boundary treatments, depending on how numerical boundary conditions are formulated and then converted together with the physical boundary conditions into a conservative form. The first assigns the boundary values or their derivatives for physical boundary conditions and computes other unknown variables (corresponding to the so-called numerical boundary conditions) mainly by space or space–time extrapolation [11], [12], [13], [10], [14], [15], [6], [8], [16], [17], [18]. In particular, for example, the subsonic outflow conditions of Thompkins and Bush [13] and Liou and Van Leer [15] fix the static pressure, and extrapolate other conservative variables. In order to rewrite the discretized boundary equation using conservative variables, nonlinear terms have to be discarded due to the high-order nonlinear terms in the energy equation. This inevitably introduces linearization errors. MacCormack [14] and Wright et al. [6] applied the extrapolate directly to the delta form of the discretized equations on the boundaries, and explicit boundary conditions are applied after each time integration. Yee et al. [10] used the characteristic projection to differentiate the physical and numerical boundary conditions, and the characteristic variables of the numerical boundary conditions are extrapolated from the interior values. While some primitive boundary values corresponding to the physical conditions are known a priori, other unknown values are computed through the reverse projection. Tidriri [8] also used the theory of characteristics in implementing far-field the inflow and outflow conditions and the impermeable surface conditions. All primitive boundary values are found by solving the locally one-dimensional characteristic relations (i.e. dropping spatial derivatives in the other two directions) and thus solutions of inner grids appear explicitly in the boundary treatments. This is also equivalent to an extrapolation of the characteristic variables or the primitive variables from the interior solutions or the far-field parameters. Lee and Lee [17] used a similar strategy but computed the defections of the conservative variables. Advancing the boundary defections along with the inner grids results in extremely complicated left implicit operators. Osusky [18] used the SBP-SAT method in a Newton–Krylov–Schur solver, where the boundary conditions are softly enforced through a simultaneous approximation term (SAT) added to the discretized NS equations using the summation-by-parts (SBP) schemes [19], [20]. In addition that the SAT has to be designed specially according to the SBP schemes, the physical boundary values are prescribed and others are also extrapolated from inner solutions to form the targeted values in the SAT. Overall, the extrapolation and various approximations impose extra constraints on the numerical boundary conditions. This may compromise the convergence rate and even stability of the implicit computations [21]. In fact, the stability is found very sensitive to the particular space–time extrapolation methods, as demonstrated by Yee et al. [10] and section 3.1 of this study. Additionally, a post-update correction is usually required after each iteration due to the inherent errors.
Without using any extrapolation for numerical boundary conditions, the second approach develops the PDE type of boundary equations, so that regular spatial and time discretizations can be applied to derive the discretized boundary equations. This essentially requires the exact governing equation for numerical boundary conditions. Except that Thompson [22] used part of the original NS equations (for example, the continuity equation for no-slip walls) for this purpose, most often a subset of the characteristic decomposed governing equations are adopted because they are otherwise challenging to design. In the characteristic decomposed form of the NS equations, Chakravarthy [23], Rai and Chaussee [24] discarded the incoming wave equations and replaced them with the time derivatives or spatial derivatives of the dependent variables that are compatible with the physical boundary conditions, for example in the form , where and p is the number of physical boundary conditions. The equations are combined with those representing outgoing characteristic waves and then transformed back to the conservative equations. Therefore, traditional one-sided spatial schemes can be used on the boundaries, and eventually advanced in the entire computational domain using the same time integration scheme along with the inner grids. Although the governing equations are expected to be accurate in terms of the physical boundary conditions, the final equations and their implicit solutions are very complicated. Due to this complexity, only very few publications adopted this strategy. Because of the errors arising from the one-dimensional approximation and inherent linearization, the boundary values are also found to drift away from the specified values and a post-update correction is thus required.
The differences and similarities between these two types of implicit boundary treatments deserve further discussion. On the one hand, the first approach usually discarded the original PDEs on the boundaries. Both the physical and numerical boundary conditions are realized in the discretized, algebraic equations where the solutions at interior grids appear. With the extra constraints introduced in the extrapolation or approximations, they may not be compatible with the original PDEs. But in the second approach, development of the boundary equations and their implicit discretizations are independent. The numerical boundary conditions are part of the original equations, to which traditional one-sided schemes may be applied. This is shown to be not only more convenient in computer code implementation, but also produce more stable and more accurate solutions as compared to the first approach. On the other hand, both approaches contain the known physical boundary values in terms of the primitive, conservative, or characteristic variables. The derived boundary equations, either in discretized or PDE forms, are complicated and remarkably different from the original PDEs. Development of some boundary conditions involving complicated combinations of the primitive variables, for instance the inlet total conditions, can be very laborious if not impossible. Boundary treatments for multidimensional problems are also very challenging and complicated when different boundary conditions are expected for the same grids at the edge or corner of the computational domain [25], [17]. Additionally, errors often arise in the linearization or approximations due to the nonlinear terms in the NS equations, this demands a post-update correction which may compromise the convergence rate of the simulations.
The purpose of this study is to introduce a new methodology for implicit boundary treatment of the NS equations. It consists of a generalized (boundary) equations and the implicit discretization techniques. It differs from the existing approaches mainly in two respects. Firstly, different from the first approach and similar to the second one, a set of boundary equations are derived in the PDE form, so that traditional one-sided spatial discretization schemes can be used directly to solve these equations along with those of inner grids. The upwind property is preserved and the numerical procedure are proved to be very stable. Secondly, with help of the characteristic decomposition but different from the second approach, it is shown that the boundary treatment actually adds a correction vector to the original NS equations. A simple correction matrix T is introduced in the NS equations, and the final PDE form is applicable for both the interior domain () and boundary conditions (). Without any approximation, the physical boundary conditions are satisfied accurately, independent of the time and spatial schemes so that any spatial discretization schemes may be used (this is different from the SBP-SAT method where the penalty-like SAT for soft boundary treatment are closely related to the SBP discretization scheme). Consequently, not only a post-update correction is unnecessary, but also the implicit solution techniques are made significantly easy to implement using, for example, the data-parallel lower–upper relation method and the Newton iteration method (combined with the GMRES method for subsidiary iterations). For multidimensional problems, different boundary conditions at the same corner or edge points can be treated successively on a direction-by-direction basis, while still showing an excellent accuracy.
It should be stressed that this study does not aim at examining the stability property of various implicit boundary treatments. At this point, no stability analysis has been made regarding these boundary equations, yet one-sided spatial schemes on the boundaries still satisfy the upwind property and thus are expected to produce stable solutions. In addition, this study does not focus on the development of high-order implicit schemes and the time accuracy is beyond the scope of this study. The numerical experiments focus only on steady flows, because they are still of great interests in application, but the proposed theory can be used directly to unsteady boundary conditions.
The rest of this paper is organized as follows. In section 2, a set of generalized equations are developed that are applicable for both the interior grids and the boundaries of a truncated computational domain, and the implicit solution techniques are described. Numerical experiments are presented in section 3 to demonstrate the accuracy and stability the proposed boundary equations. Finally, conclusions and future work are discussed in section 4.
Section snippets
Implicit boundary equations for the NS equations
As discussed in section 1, one may apply space–time extrapolation or one-sided spatial schemes on the boundaries for numerical boundary conditions. However, the former may degrade the computational stability and accuracy, and the latter requires specially designed boundary equations because the original NS equations do not satisfy all boundary conditions automatically. Both approaches yield complicated discretized boundary equations and are challenging to implement in computer code.
This section
Numerical results
In this section, four benchmark cases are simulated for validation purpose. The first solves the nonlinear Burger's equation with a steady state solution, in order to demonstrate the advantage of direct time–space discretization over space–time extrapolation in terms of the computational accuracy and stability for implicit boundary treatment. The second and the third cases demonstrate the accuracy and stability of the proposed boundary treatment for NS equations, and the last shows its improved
Conclusions
Implicit boundary treatment is a particularly difficult problem in modern CFD simulations. Special treatments are imposed on the boundaries due to the absence of flow solutions beyond the boundaries and the yet unknown solutions at future times. Different boundary treatment may not only affect the computational accuracy, but also significantly alter the eigenvalues, and consequently the stability, of the entire system. In general, the implicit boundary treatment should be formulated in the
Acknowledgments
The author would like to thank Dr. John A. Ekaterinaris from Embry-Riddle Aeronautical University, USA for many helpful discussions. Financial support from the National Natural Science Foundation of China (grant number 11772262), Fundamental Research Funds for the Central Universities (grant number 3102018zy002) and the 111 project of the Department of Education of China (grant number B17037) are greatly appreciated. Anonymous reviewers are also acknowledged for their constructive suggestions
References (38)
Preconditioning techniques for the Newton–Krylov solution of compressible flows
J. Comput. Phys.
(1997)- et al.
Jacobian-free Newton Krylov methods: a survey of approaches application
J. Comput. Phys.
(2004) - et al.
On the application of boundary conditions to time dependent computations for quasi one-dimensional fluid flows
Comput. Fluids
(1977) - et al.
Boundary treatments for implicit solutions to Euler and Navier–Stokes equations
J. Comput. Phys.
(1982) - et al.
Time-stable boundary conditions for finite difference schemes for solving hyperbolic systems: methodology and application to high-order compact schemes
J. Comput. Phys.
(1994) - et al.
Review of summation-by-parts schemes for initial–boundary-value problems
J. Comput. Phys.
(2014) - et al.
Stability analysis of numerical boundary conditions and implicit difference approximations for hyperbolic equations
J. Comput. Phys.
(1982) - et al.
Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows
J. Comput. Phys.
(2008) Time-dependent boundary conditions for hyperbolic systems, II
J. Comput. Phys.
(1990)- et al.
Boundary conditions for direct simulations of compressible viscous flows
J. Comput. Phys.
(1992)