Abstract

A computational fluid dynamics (CFDs) method utilizing unstructured grid technology has been employed to compute vortical flow around a delta wing with sharp leading edge, which is specially known as the geometry of the second international vortex flow experiment (VFE-2). In VFE-2, delta wings with different leading edges had been broadly investigated by experiments, which resulted in a special database for CFDs codes validation. The emphasis of this paper is to investigate the effectiveness of an adjoint-base grid adaptation method for unstructured grid in capturing concentrated vortices generated at sharp edges or flow separation lines of lifting surfaces flying at high angles of attack. Earlier experiences in vortical flow simulations had indicated that the vortex behavior is highly dependent on the local grid resolution both on body surface and space. The adjoint-based adaptation method used here is hoped to save grid points with a reasonable grid resolution for vortical flow simulations. The basic idea is to construct a new adaptive sensor in a grid adaptation process with the intent to tell where the elements should be smaller or larger by introducing an adjoint formulation to relate the estimated functional error to local residual errors of both the primal and adjoint solutions.

1. Introduction

In recent years, computational fluid dynamics (CFDs) has matured as a routine tool to predict complex genuine flows at designed condition. However, it is still difficult to deal with complex flows at off-design conditions where flow separations and vortices are dominated. For these flows, inappropriate grid resolution has been recognized as one of the dominant reasons for inaccurate predictions of flow features and aerodynamic forces. Actually, in general, the computational grid may become rapidly coarser as it becomes far away from the body surface. In this case, the vortices would be highly diffused due to the numerical discretization error. A possible way to minimize the numerical dissipation of the vortices is to use a highly dense grid. But flow computations around three-dimensional complex bodies with large-scale separations and vortices are difficult within the realistic number of grid points. In order to obtain numerical results for such flows with a given level of accuracy by a limited capacity of the computer, the method for grid adaptation is very useful for practical computations, which offers the possibilities to avoid the use of overly refined grids to guarantee accuracy. The basic premise is to locally enrich the computational grids in regions which most adversely affect the accuracy of the final solution. Its emphasis is to define a proper adaptive criterion (also known as the adaptive sensor) combined with distinct methods to automatically generate an adaptive mesh.

Currently, many sensors for unstructured grid adaptation are available, for example, the wake sensor based on difference in velocity, the shock sensor based on ratio of total pressure, the vortex sensor based on creation of entropy, and the vortex sensor based on eigenvalue analysis of the velocity-gradient tensor, and so forth [1]. These adaptation methods are often referred to as feature-based adaptation, which mainly focus on resolving discontinuities or strong gradients in the flow field. The assumption herein is that regions of larger gradients are associated with regions of larger error. The use of unstructured grid technology for CFDs simulations allows more freedom in adapting the discretization of the meshes to improve the fidelity of the simulation. Unfortunately, continuous local refinement of the dominant features of the flow does not necessarily guarantee that the global error will simultaneously be reduced. In some cases, this may even lead to incorrect results [2]. The reason is that the adaptive sensor has failed to correctly indicate the location where the elements should be larger or smaller due to errors elsewhere in the flow field. Also, since the feature-based adaptation methods are nearly without regard to an output quality (e.g., lift coefficient, drag coefficient, moment coefficient, and so on), they may distort to resolve the flow in regions with large local error that may have a minimal effect on the output.

In this paper, an adjoint-based adaptation with isotropic h-refinement for unstructured grid is introduced to catch the vortex-dominant flow. These kinds of methods are firstly proposed in finite element communities [3] and were then extended to finite volume computations, for example, [410]. The basic idea of the adjoint-based adaptation method is to construct an adaptive sensor from a more robust and accurate error estimation by introducing the adjoint formulation to relate the estimated functional error to local residual errors of both the primal and adjoint solutions. The adjoint solutions provide a very powerful approach to compute output error estimation as well as to systematically adapt grids to reduce spacial discretization errors. It has also been observed that the adjoint-based adaptation can greatly reduce unnecessary grid point compared with the feature-based adaptation both for 2D and 3D flows [11]. The application of the adjoint-based adaption here is hoped to economize the grid points with a reasonable grid resolution for vortical flow simulations. The test case is a swept delta wing with sharp-leading edge from the second interactional vortex experiment (VFE-2) whose experimental data have been authorized to use by German Aerospace Center (DLR), Göttingen, Germany. This test case is also from NASA Langley national transonic facility (NTF) that has been used to illustrate effectiveness of the vortex sensors developed by us [1]. Besides, inviscid transonic flows over RAE2822 airfoil and ONERA M6 wing are given to validate the adjoint-based adaptation method used here.

2. Methodology

2.1. Flow Solver

The flow solver used in this paper is WoF90, which is developed by ACTRI. It is a solver employing hybrid unstructured grid both for Euler and Reynolds-averaged Navier-Stokes (RANS) equations. The solver can deal with arbitrary types of unstructured elements including tetrahedrons, hexahedrons, prisms, and pyramids. It is based on an edge-based formulation and uses a node-centered finite volume technique to solve the governing equations. The control volumes are nonoverlapping and are formed by a dual grid, which is computed from the control surfaces for each edge of the primary input mesh. The equations are integrated explicitly towards steady state with Runge-Kutta time integration. The spatial discretization is either central with artificial dissipation or upwind; both approaches are second order accurate. The solver also adopts an agglomeration multigrid algorithm and an implicit residual smoothing algorithm to accelerate the convergence of a simulation.

In this current study, Euler equation is solved by WoF90 for conserved variables , where is density, is velocity components, is total energy per unit volume. The semidiscretization of Euler equation in space using finite volume technique leads to the following ordinary differential equation (ODE): where is the flow equation residual to steady-state solution and is volume of a control volume.

2.2. Adjoint Solver

Let is the objective output function we are interested in. Often we can choose lift coefficient, drag coefficient, moment coefficient, or their combinations, and so on, as the objective output quality. The adjoint variable is defined as the linearized effect of a flow residual source term on the output function [10, 11]: It can be found by solving a set of linear equations, which is derived with application of the chain rule for differential After the flow solution is known, this set of linear equations can be solved much like the flow solver by adding a pseudotime item in order to result in a more robust adjoint solver [7, 12] where is the adjoint equation residual This equation has a form similar to the flow equation so that it can be integrated toward steady solution with the same Rung-Kutta time integration as the solver. An inviscid adjoint solver has been implemented in WoF90.

2.3. Adaptive Sensor Based on Error Estimation with Adjoint Solution

The adjoint solutions can provide an efficient way to estimate error of the objective output function . Following the derivation in [79], we give a brief derivation here. Firstly, let us assume is a solution on the initial coarser mesh, which can be defined as a perturbation to a solution on a much finer mesh as where the fine-mesh solution and the error are not known in practice. The resulting errors in the output function and the flow equation residual are, respectively, Here, the residual on the much finer mesh is assumed to be zero, that is, . Linearizing about the fine-mesh solution yields, From (8) and by introducing the adjoint variable from (3), we can give an alternative formula for the error in the output function as where is the adjoint solution on the much finer mesh, which is often undesirable. Instead, it is assumed an approximate adjoint solution is available as a substitute. In a notation of the error of adjoint solution , equation (10) can now be rewritten as where is the adjoint equation residual on the much finer mesh, and is assumed.

As shown in (11), the error in the output function can be classified two classes: a computable correction that can be evaluated with the solution on current coarse mesh and a remaining error that generally cannot be evaluated without solving for qualities on the much finer mesh, that is, The average of the absolute values of the two remaining errors is defined as the adjoint-based adaptive sensor where, in practice, both the adjoint solution error and the flow solution error would be replaced by interpolation errors and are obtained by interpolating into a globally refined mesh.

2.4. Grid Adaptation and Smoothing

A general introduction to the flow chart of grid adaptation can be found in [1, 13]. In the current study, the local h-refinement technique is used due to its robustness and efficiency, which divides existing elements in the mesh according to initial solution (also see examples in [14, 15]). For an inviscid calculation, the unstructured grid only has elements of the tetrahedrons in three dimension (3D) or of the triangles in two dimension (2D). These elements are refined isotropically employing the binary division according to marked edges. It should be pointed out that since the binary division will bring the so-called “hanging node" and unavoidable elements differ from tetrahedrons and triangles, the method to eliminate these nodes or elements is also used here [16]. Besides, after adaptation the mesh may have a worse quality (e.g., the aspect ratio) for some of elements, while they may significantly impact the accuracy and convergence of a solution. We recommend to import the adapted mesh into a general grid generator (e.g., ICEM-CFDs, etc.) for smoothing the mesh in order to improve the global mesh quality.

To use the adjoint-based adaptive sensor defined in (13) and (14) to indicate where to be refined, a globally h-refinement mesh is firstly created by Octree division, in which a new node is inserted at the midpoint of each existing edge, and then both solutions of the flow equation and the adjoint solution on the much finer mesh are evaluate by piecewise linear interpolation. The local adjoint-based adaptive sensor now can be evaluated from (13) by exploring all nodes of the much finer mesh within each element The criterion for refinement is given as follows: given a user-specified desired error , let , where is the total number of elements in current mesh, and judge if is satisfied. In this paper, a desired error is used for all cases.

3. Results

3.1. Validation for 2D Case

Inviscid transonic flow over RAE2822 airfoil is chosen to validate the rationality of adjoint-based adaptive method for two-dimensional test cases. This test case is a classical validation case, which has been used worldwide, especially in EUROVAL project and by AGARD. The simulation is performed at the flow condition as and , which is a corrected condition from case 9 in EUROVAL project.

The initial mesh has a total of 15,852 nodes and 31,388 triangles, which is shown in Figure 1. For flow solution, the central scheme is employed for spatial discretization, and 3 multigrid levels with W-cycle are used for acceleration. The result of Mach contours employing the initial mesh is given in Figure 2. The primary flow feature is evident.

To apply adjoint-based adaptation method, we choose the lift coefficient as the objective output function for adjoint solutions. The final adapted mesh using adjoint-based sensor and isotropic adaptation with h-refinement is shown in Figure 3, which contains a total of 20,763 nodes and 41,141 triangles after mesh adaptation. It is clearly shown that a large portion of new nodes have been appended to leading edge, upper surface and shock regions. This finding also fits the primary feature of the flow very well. The initial mesh has been sufficiently refined in high-difference regions in flow velocity.

3.2. Validation for 3D Case

Inviscid transonic flow over ONERA M6 wing is chosen to validate the rationality of adjoint-based adaptive method for three-dimensional test cases. This configuration is simulated at the flow condition at and . The initial mesh contains a total of 52,144 nodes and 296,517 tetrahedrons, where there are 5,459 boundary nodes on the wing surface. The initial wing surface mesh is shown in Figure 4, which has a coarse spacing, especially at the trailing edge. It is designed with the intention to resolve the surface curvature of the leading edge and the wingtip, but it has no particular considerations of flow features. For the flow solution, the central scheme is used for spatial discretization, and 3 multigrid levels with W-cycle are used for acceleration.

Again, we choose the lift coefficient as the objective output function for adjoint solutions. After grid adaptation using adjoint-based sensor and isotropic adaptation with h-refinement, the adapted mesh contains a total of 128,386 nodes and 740,366 tetrahedrons, where there are 8,041 boundary nodes on the wing surface. The adapted wing surface mesh is shown in Figure 5. A dense region with distinct lambda structure is found in the mesh after adaptation, which just indicates the shock region of M6 wing under current circumstance. It is also found that sufficient nodes have been enriched both at the leading edge and the trailing edge.

The Mach contours on the adapted upper wing surface are shown in Figure 6. A lambda shock structure is clearly visible in these Mach contours, which also indicates the flow physics of M6 wing under the current flow condition. Figure 7 further gives the comparisons of the surface pressure coefficient (noted by ) at different typical location along the spanwise direction. The experimental data are compared with the original results and the adapted results for both upper and lower surface of the wing. It is evident that the surface pressure distribution has been successfully improved by the grid refinements with current adjoint-based adaptation method.

3.3. Application to Complex Vortical Flows of Delta Wing

As mentioned above, the test case of delta wing with sharp leading edge is from VFE-2 as well as from NASA NTF. The emphasis of this test is to validate the vortex-capturing capability of the adjoint-based adaptation method. It is worth mentioning that while it is hard to accurately predict this kind of complex vortex flows by Euler calculation, the knowledge of vortex location obtained with grid adaptation from the inviscid solution may greatly help viscous simulations [17].

The detailed description of the model can be found in [18]. The experiment was conducted both in the NASA NTF [18] and in the transonic wind tunnel at Göttingen (TWG) of German Dutch wind tunnels (DNW) [19]. In the current study, all computations were performed on a semispan model at the flow condition at and . For the flow equation solution, the central scheme is used for spatial discretization and 3 multigrid levels with W-cycle is used for acceleration.

The initial mesh contains a total of 64,835 nodes and 365,681 tetrahedrons, where there are 9,464 boundary nodes on the wing surface. No attempt has been made to cluster extra grid points in the field at the vortex location. The initial surface mesh on the upper surface is shown in the upper portion of Figure 8. A fine grid resolution has been prescribed at the edges of the wing, whereas the mesh on the flat portion of the wing is nearly uniform. The corresponding flow solution on this mesh is displayed in the lower portion of Figure 8, where the distribution of pressure coefficient is represented by variation of colors. The footprint of a leading-edge vortex on the surface is clearly found from the low-pressure region.

To apply the adjoint-based adaptation, we also choose the lift coefficient as the output function. After the application of grid adaptation, the adapted mesh contains a total of 139,682 nodes and 766,948 tetrahedrons, where there are 13,333 nodes boundary nodes on the wing surface. The adapted upper surface mesh and the corresponding Euler solution on the adapted mesh are shown in Figure 9. It is evident that a large portion of new nodes have been enriched on the leading edge, the trailing edge, and the concentrated primary vortex regions such that the footprint of a leading-edge vortex on the surface becomes clearer. Furthermore, the primary vortex appears to be shifted slightly towards the wing root as compared with unadapted solution. In Figure 10, two typical sections at and were chosen for grid resolution comparison before and after mesh adaptation. Obviously, the mesh has been refined in the vortex region above the upper wing surface, while remaining unchanged under the lower wing surface. It can be also noticed that the grid density of refined regions varies with movement of the vortex core. However, it is also been found that, unlike using the adaption based on some elaborated vortex sensors [1], the spatially refined elements from adjoint-based adaptation were not over concentrated along the vortex core.

In Figure 11, the spanwise surface pressure coefficient distributions at different chordwise positions (, and 0.8, resp.) were compared between before and after mesh adaptation. In these figures, experimental data both from NASA and DLR are depicted. The improvement for prediction of the vortex flows is obvious by using the adjoint-based adaptation and h-refinement. The unadapted Euler simulation has produced a poor solution with underpredicted vortex in the upstream stations due to uniform grid spacing, especially, it has caused a severe shift of the peak value of on the upper surface to the wing tip. After adaptation, things have been changed greatly. At all chordwise positions except , it appears to give a more reasonable result. At position, similarly poor results are obtained with initial and adapted meshes. It may be said that the two meshes have the same insufficient grid resolution around this section (compare Figure 8 and Figure 9). This deficiency may be ameliorated by adapting the new mesh or by altering the desired error level used in the adjoint-based adaptation.

4. Conclusions

An adjoint-based adaptation on unstructured grid has been implemented and is applied to simulate the complex vortex flow of a delta wing with sharp leading edge at high angle of attack. Firstly, an adaptive sensor for grid adaptation is constructed from an accurate error estimation to output function with adjoint solution, which combines residuals of the flow equation and the adjoint equation. Then, the rationality and capability of the adjoint-based adaptation is validated by choosing the transonic flow over RAE2822 airfoil as 2D case and the transonic flow over ONERA M6 wing as 3D case. Reasonable results have been obtained. After that, results from initial numerical investigation of delta wing are given. We can conclude that although the spatial refinements from adjoint-based adaptation were not over concentrated along the vortex core of the flow, the results were greatly improved. The greatest advantage of an adjoint-based adaptation is that it easily reaches an improved and convergent grid resolution and saves unnecessary nodes which have a minimal effect on the desired output function.

It should be pointing out that although all simulations in this study are limited to inviscid Euler computations, the adjoint-based adaptive method itself can be extended to viscous calculations by developing a viscous adjoint solver and using anisotropic refinement for grid within boundary layers [9]. These work is currently under way at ACTRI.

Acknowledgments

This work is supported by the State Key Development Program of Basic Research of China (973) under Grant no. 2009CB723804. The authors would like to thank the German Aerospace center (DLR) for proving and authorizing to use the experimental data of delta wing in VFE-2 project. Mr. YANG Zhenhu, one of our colleagues at ACTRI, is also appreciated due to his part of work on adjoint-based mesh adaptive method.