Elsevier

Journal of Computational Physics

Volume 375, 15 December 2018, Pages 641-658
Journal of Computational Physics

Implicit boundary equations for conservative Navier–Stokes equations

https://doi.org/10.1016/j.jcp.2018.09.010Get rights and content

Highlights

  • An implicit boundary treatment consisting of a generalized boundary equation in PDE form and its implicit solution techniques is proposed.

  • A correction matrix T is introduced in the original NS equations to realize the desired boundary conditions.

  • One-sided spatial schemes can be directly applied to discretize the conservative equations on boundary points, and are still shown stable.

  • The boundary residuals are computed more accurately, and accelerate the convergence rates of implicit solutions significantly.

  • Computer code implementation is made significantly easy, especially for multidimensional problems.

Abstract

Co-existence of the physical and numerical boundary conditions makes implicit boundary treatment a particularly difficult problem in modern CFD simulations. Previous studies adopted space–time extrapolation or specially designed partial differential equations on the boundaries that are different from those of interior points. They are often formulated in terms of primitive variables, and are very challenging for complicated boundary types to be converted to, and implemented in, the conservative variables that are preferred in numerical simulations. More importantly, different boundary equations or different extrapolation techniques may compromise the stability, accuracy, or convergence rate of the A-stable schemes that are developed for interior points. A new methodology for implicit boundary treatment is proposed in this study. By introducing a simple correction matrix T, a set of generalized equations Q/t=(I+T)R are developed in terms of conservative variables. It is applicable for both the interior domain (T=0) and the boundaries (T0). It is in a partial differential equation form, satisfies the boundary conditions accurately, independent of the time and spatial discretizations. Any one-sided schemes can be used on the boundaries but still maintain the upwind property. Implicit solution techniques are made significantly easy to implement using, for example, the data-parallel lower–upper relation method and the Newton method (combined with the GMRES method for subsidiary iterations). Numerical experiments show that the proposed methodology produces stable simulations for very large CFL numbers and preserve the imposed boundary values accurately.

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:Qt=R(Q)=EˆξFˆηGˆζ+Eˆvξ+Fˆvη+Gˆvζ 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:ut=R(u)=fx where, R is again referred to as the residuals involving spatial derivatives but in the x direction only, and the flux is f(u)=u2/2. The numerical solutions are sought at discretized time instances tn=nΔt (n=1,2,3...) and grid points xi (i=1,2,...,N) within axb, 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:uin+1uinΔt=Rin+1=fx|in+1

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 tn using the Taylor series expansion, and organize the resulting equation into a delta form:(1ΔtRu|n)Δun=Rn where, Δuin=uin+1uin and the numerical solutions will be uin+1=uin+Δuin. Or alternatively, one may re-write Eq. (3) into F(qi)=(qiuin)/ΔtR(qi)=0, and then solve it using the Newton iteration:(1ΔtRq|ik)Δqik=[qikuinΔtRik]qk+1i=qik+Δqik where, the iteration can start with qi0=uin, and accurate solutions are recovered as uin+1=qik when F(qk) 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:(1ΔtRu|n)Δun=RHS(tn,un) where, the left implicit operator (in particular, the Jacobian matrix R/u) depends on the spatial discretization schemes. For example, when finite difference schemes are used,Rin+1=fx|n+1i=j=lrci+jfi+jn+1 where, l and r indicate the left and right shifts of the grid points in the stencil. Therefore,Ru|inΔuni=j=lrci+jfu|i+jnΔui+jn and with f/u=u, 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 u>0 in the entire domain axb, u=ua at the inlet x1=a must be given to establish the problem, and this is referred to as the physical boundary conditions. At the outlet xN=b, 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, uN=uN1, 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 B(Q)=0 are preferred, so that they can be discretized intoBQ|nΔQn=B(Qn) 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 Bi(Q)=0, where i[1,p] 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 Q/t=(I+T)R is applicable for both the interior domain (T=0) and boundary conditions (T0). 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)

  • C.K.W. Tam et al.

    Dispersion-relation-preserving finite difference schemes for computational aeroacoustics

    J. Comput. Phys.

    (1993)
  • C. Hirsch

    Numerical Computation of Internal and External Flows

    (1990)
  • J. Nordstrom et al.

    Well-posed boundary conditions for the Navier–Stokes equations

    SIAM J. Numer. Anal.

    (2005)
  • S. Yoon et al.

    Lower-upper implicit schemes with multiple grids for the Euler equations

    AIAA J.

    (1987)
  • S. Yoon et al.

    A lower–upper symmetric-Gaussian Seidel method for the Euler and Navier–Stokes equations

    AIAA J.

    (1988)
  • M.J. Wright et al.

    A data-parallel LU relaxation method for the Navier–Stokes equations

  • M.J. Wright et al.

    Data-parallel line relaxation method for the Navier–Stokes equations

    AIAA J.

    (1998)
  • M.D. Tidriri

    Krylov Method for Compressible Flows

    (1995)
  • H.C. Yee et al.

    Boundary approximations for implicit schemes for one-dimensional inviscid equations of gasdynamics

    AIAA J.

    (1982)
  • Cited by (0)

    View full text