A fast and automatic full-potential finite volume solver on Cartesian grids for unconventional configurations

Fanxi LYU,Tianhang XIAO,Xiongqing YU

College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China

A fast and automatic full-potential finite volume solver on Cartesian grids for unconventional configurations

Fanxi LYU,Tianhang XIAO,Xiongqing YU*

College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China

Available online 20 April 2017

*Corresponding author.

E-mail addresses:lvfanxi@nuaa.edu.cn(F.LYU),xthang@nuaa.edu.cn(T.XIAO),yxq@nuaa.edu.cn(X.YU).

Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

http://dx.doi.org/10.1016/j.cja.2017.03.001

1000-9361©2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

To meet the requirements of fast and automatic computation of subsonic and transonic aerodynamics in aircraft conceptual design,a novel finite volume solver for full potential flows on adaptive Cartesian grids is developed in this paper.Cartesian grids with geometric adaptation are firstly generated automatically with boundary cells processed by cell-cutting and cell-merging algorithms.The nonlinear full potential equation is discretized by a finite volume scheme on these Cartesian grids and iteratively solved in an implicit fashion with a generalized minimum residual(GMRES)algorithm.During computation,solution-based mesh adaptation is also applied so as to capture flow features more accurately.An improved ghost-cell method is proposed to implement the non-penetration wall boundary condition where the velocity-potential of a ghost cell is modified by an analytic method instead.According to the characteristics of the Cartesian grids,the Kutta condition is applied by specially computing the gradients on Kutta-faces without directly assigning the potential jump to cells adjacent wake faces,which can significantly improve the solution converging speed.The feasibility and accuracy of the proposed method are validated by several typical cases of sub/transonic flows around an ONERA M6 wing,a DLR-F4 wing-body,and an unconventional figuration of a blended wing body(BWB).The validation cases demonstrate a fast convergence with fully automatic grid treatment and computation,and the results suggest its capacity in application for aircraft conceptual design.

©2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Cartesian grids;

CUT-cell;

Finite volume method;

Full potential equation;

Grid adaptation;

Kutta condition;

Non-penetration condition

1.Introduction

In aircraft conceptual design and multidisciplinary optimization,numerical methods for aerodynamic analysis are replacing traditional engineering methods1,2due to the fact that numerical methods have higher accuracy and can be applied to unconventional aircraft configuration studies.Generally,various design concepts are required to be evaluated in the conceptual design stage,and each design concept may need a multidisciplinary optimization,thus resulting in a large number of calls for aerodynamic analysis;therefore,a numerical method used in this stage must be rapid and fully automatic to meet the requirements.

In terms of accuracy and efficiency,full potential solvers,among various numerical methods,are regarded as the most suitable ones for subsonic and transonic flow analysis in the conceptual design stage.Though methods based on Euler equations and Navier-Stokes equations have high accuracy and versatility,the huge computational cost makes them difficult to apply in conceptual design.Panel methods can respond to aerodynamic analysis calls rapidly,1but cannot calculate nonlinear flows present at transonic speeds and are limited to either subsonic or supersonic flow fields.Full potential methods are capable of accurately computing subsonic and transonic flow fields at a relatively low computational cost,and therefore,they can meet the requirement of rapid response.1,3Three methodologies are documented that can solve full potential flows,i.e.,finite difference method(FDM)on structured grids,finite element method(FEM)on unstructured/Cartesian grids,and finite volume method(FVM)on structured/unstructured grids.

The FDM for solving the full potential equation was developed in 1970s and has reached a certain level of maturity,e.g.,FLO 22,4one of the most successful tools of this type,has been widely used for years for wing design at Boeing Company.5However,it is only on structured grids that the FDM can be applied to solve the potential equation,but often structured grids are very difficult to generate;therefore,the applications of the FDM are limited to single wing or simple wing-body configurations.In addition,automatic implementation of this method is a very difficult task,particularly for complex geometries,due to the barriers of structured mesh generation.

The FEM subdivides a flow field into a number of small finite elements,and the governing equation of the potential flow is modeled on these elements to approximate a solution by minimizing an associated error function.An attractive advantage of the FEM is its ability to handle complicated domains and arbitrary shapes of elements,and therefore it can be applied on either unstructured or Cartesian grids.The most representative examples of this method is TranAir,6,7a full potential solver,which employs the FEM to solve the governing equation on Cartesian grids.The unique features of Cartesian grids,e.g.,strong capability of space-filling,automatic mesh generation,and easy solution adaption,make TranAir a robust tool for conceptual design,and it has seen heavy use in industrial applications.However,one problem related to the FEM is that it requires special care to ensure a conservative solution,particularly for flow fields with discontinuities,e.g.,transonic flows with shock waves.Additionally,the FEM may consume more memories and have slower solution times than the FDM and the FVM due to its specialties in numeric aspect.Regarding automatic implementation,the capability of the method is probably restricted by the need of additional treatments on geometries or boundaries,e.g.,TranAir has to run an additional process module,aero grid and paneling system5(AGPS),to represent an input geometry into a set of quadrilateral panels for the convenience of subsequent computation.5

During the past decades,a few computational tools have been developed based on the FVM to solve a nonlinear full potential equation.In the FVM,the potential equation is recast in a conservative form and then solved over discrete control volumes.A natural advantage of this discretization is that solution conservation can be guaranteed in principle.Furthermore,the FVM benefits from its less memory usage and faster solution speed.Currently,the FVM for potential flows is implemented mainly on structured or unstructured grids.Jameson and Caughey8–10and Bolsunovsky et al.11contributed great efforts on the FVM on structured grids,and the full potential solvers,FLO 27/28/30 and BLWF developed by them respectively,have been successfully applied for aerodynamic analysis of wings or wing-body figurations in engineering.However,the problems related to structured grids,as mentioned before,make them rapidly inapplicable for unconventional configurations.Later on,attempts were also made by Neel12to apply the FVM on unstructured grids where a solution for a complicated geometry thus becomes possible.Due to its capability of space filling,an unstructured solver is superior to a structured solver in terms of robustness for unconventional configurations.However,it is still impractical for an unstructured solver to fully automatically run a solution since the automatic generation of an unstructured surface mesh is still a difficult task and the generation of a volume mesh is also at a relatively high level of computational cost.

Compared to structured or unstructured grids,Cartesian grids have a strong capability of space-filling and are suitable for complex geometries,and moreover,it can be implemented automatically with less difficulty for mesh generation and solution adaptation.13,14These features make Cartesian grids an ideal way for rapid and automatic aerodynamic analysis for various geometric changes in aircraft conceptual design.With respect to spatial discretization,the FVM provides a reliable choice to ensure a conservative,accurate,and efficient solution because of its distinguishing features.Consequently,it is wellreasoned to combine the FVM with the Cartesian grid technique for automatic full potential analysis at each design stage.However,in spite of the significant progress in computational techniques for full potential flows,methodologies related to the FVM on Cartesian grids have not yet been properly addressed,especially when both commonality for unconventional designs and fully automatic process are strongly stressed.

In the present paper,a full potential finite volume solver on Cartesian grids is developed for the sake of providing a new fast and automatic transonic aerodynamic analytical tool in the conceptual design stage particularly for unconventional configurations.To our knowledge,this work is the first attempt to implement the FVM on Cartesian grids for full potential flows.A few innovative strategies and methods related to the discretization scheme,wall boundary condition,and Kutta condition are proposed in the circumstance that techniques used on structured or unstructured grids may suffer poor performance when they are migrated on Cartesian grids.The developed solver and proposed methods are validated by several typical cases including complicated unconventional figurations.The results demonstrate a fast convergence history with highly automatic grid treatment and solution process,suggesting its capacity in application for aircraft conceptual design.

The remainder of the paper is organized as follows:Section 2 brief l y introduces the techniques of automatic Cartesian grid generation.The finite volume discretization and implicit solution for full potential equations are presented in Section 3.Boundary conditions and their corresponding implementation,including the wall boundary,Kutta condition,and far-field condition,are given in Section 4,followed by validations and results in Section 5 and some conclusions in Section 6.

2.Cartesian grid method

This section describes the automatic generation of Cartesian grids.The Cartesian grid method is implemented in the following four main steps:(A)importing a CAD model and organizing geometric information by an efficient data structure;(B)applying geometric refinement on the initial Cartesian grids;(C)classifying and smoothing the grids;and(D)applying the cell-cutting method to process the surface boundary with a cell-merging algorithm on the grids.

2.1.Geometric model and data structure

The input geometric model is imported by a STereoLithography(STL)format which is one of the most popular CAD formats and supported by almost any CAD tools.The STL format describes an unstructured triangulated surface by the unit normal vector and vertices(ordered by the right-hand rule)using a three-dimensional Cartesian coordinate system.These triangular facets are subsequently organized by an alternative digital tree(ADT)15data structure so that the geometric search process can be performed efficiently.

The present study applies an octree data structure to determine gird connectivity.A Cartesian grid is generated from a root cell with a depth of 0 in the octree,which represents the complete computational domain.The grid cells in the octree are recursively subdivided according to the size requirement.In the octree,each cell includes the following information:parent’s address,eight children’s addresses,cell center,cell size,octree depth,and type,as illustrated in Fig.1.Although the data are organized by the octree structure,all the grid cells are stored in successive memory.

2.2.Grid initialization and geometric adaptation

The initial uniform grids,which are generated from the root cell of the octree and refined recursively until the specif i ed global size is reached,are not fine enough to describe the wall boundary accurately.Therefore,a geometric adaptation13,14is then implemented to refine the cells near the wall boundary.A cell intersecting with the wall boundary will be refined if it meets any of the following conditions:(A)its depth is less than its minimum depth;(B)its size,the cubic root of its volume,is larger than the minimum of the sizes of the triangular facets which intersect with it;(C)within it,the angle between the normal vectors of the intersected triangular facets is larger than a certain threshold;or(D)in two adjacent cells,the angle between the normal vectors of the intersected triangular facets is larger than a certain threshold.Taking a 2D problem as an example,Fig.2 illustrates the Cartesian grids around an NACA0012 airfoil after a geometric adaptation has been performed.

2.3.Cell classifying and smoothing

To meet the requirements of meshing and flow solution,Cartesian grid cells are classified into three categories:SOLID-cells that wholly lie within solid walls,FLUID-cells being fully submerged in fluid,and CUT-cells which intersect with the wall boundary.All the CUT-cells are firstly sorted out by using an ADT data search algorithm,and then the SOLID-cells and FLUID-cells can be identified rapidly by a painting algorithm14once any one of the SOLID-cells or FLUID-cells is distinguished inside or outside the wall boundary by applying an inside-outside check technique,e.g.,the ray-casting algorithm.14

After cell classifying,a smoothing procedure is implemented to refine some cells,and therefore a much smoother variation of cell sizes can be obtained.This procedure is beneficial to improve mesh quality and simplify the follow-up process.In general,the smoothing procedure is applied to cells that meet any of the four conditions:

(1)The tree depth differs by larger than 1 between adjacent cells.

(2)A cell has children,but all its neighbors in the same tree depth do not have children.

(3)A cell does not have children,but all its neighbors in the same tree depth have children.

(4)The tree depth of a CUT-cell is larger than the depth of its FLUID-cell neighbors.

It should be mentioned that the SOLID-cells identified in this step are excluded from computation,and the FLUID-cells and CUT-cells need to be further processed in the following procedure.

2.4.Cutting cells

As the generated Cartesian grids do not conform to the wall boundary,additional treatment on the CUT-cells is needed so that the wall boundary condition can be properly implemented.Due to its accuracy of boundary and degree of development,the cutting method13,14,16,17is employed in the present study to cut the CUT-cells into boundary cells.In general,a wall boundary cell is formed by three kinds of faces:interior faces that do not intersect with the wall boundary,the wall faces,and the interior faces which intersect with the wall boundary.The first kind can inherit directly from the CUT-cells by the tree structure,the second kind is produced by reentrant polygon clipping18from the wall geometries,and the third kind is generated from the information of the tree structure and the edges which have been generated.

In terms of the treatment of the CUT-cells,it is not convenient for the octree structure to store the topological relation of the CUT-cells due to the fact that the intersection between the CUT-cells and the wall boundary is much complicated as illustrated in Fig.3.For example,in Fig.3(a),a CUT-cell may be divided into multiple volumes.The present study employs an additional ‘cell-to-face” unstructured data structure16,19to store the topology of the processed grids,where pointers are created from cell interfaces to their adjacent cells as the basic connectivity information.Therefore,different kinds of complex intersection circumstances could be processed and topological information could be stored correctly,as seen in Fig.3(c).During the above process,some boundary cells with very small volumes may be generated.To avoid any instability problem caused by these small cells,a cell-merging algorithm17is applied to merge a small cell with its largest neighbor,as illustrated in Fig.3(b).

3.Full potential finite volume solver

3.1.Governing equations

The three-dimensional steady full potential equation is given as

whereVis the velocity vector and also the gradient of the velocity potential φ,u,vandware the components of the velocity of the fluid,ρ is the density which can be computed via the isentropic formula as

where γ is the ratio of specific heats andMa∞is the freestream Mach number.

3.2.Finite volume formulation

3.2.1.Spatial discretization

The integral form of the governing equation derived by the divergence theorem is written as

where Ω is the control volume,bounded by the closed surface∂Ω,andSis the face of this control volume.By using the finite volume scheme,the integral equation is spatially discretized on each cellias

where ρjandVjare the density and velocity at facejwhich belongs to celli,respectively.njand ΔAjare the unit normal vector and the area of facej,while Resiis the residual of celli.

3.2.2.Flux computation

(1)Potential gradient reconstruction

The potential gradient,viz.,Vjin Eq.(5),is reconstructed by a K-exact least squares method.12As illustrated in Fig.4,for a given interior interface with two adjacent grid cells,all the face-neighbor cells of these two adjacent cells are chosen as a stencil for the reconstruction of the potential gradient on this interface.The potential in each cell is assumed linear as

Then the error associated with each cell can be written as

The sum of the squares of the weighted error is obtained as

whereriis the weighting factor computed from the inverse of the distance between the cell center and the face center,andmis the number of cells in the stencil.According to the least squares method,the equations are derived from the partial derivatives ofFwith respect to the unknowns as

Therefore,the solutions of Eq.(9)is written as

whereis the element of the inverse matrix of the coefficients matrix in Eq.(9),andC1,0,0,C0,1,0,andC0,0,1are the potential gradient,viz.,the velocity componentsu,v,andw,respectively.

(2)Artificial viscosity

For the consideration of numerical stability,a density biasing method is applied following the work of Hafez et al.20,where an artificial viscosity is added to the equation with the density in Eq.(1)being replaced by~ρ as

where Δx,Δy,and Δzare twice the distance from the face center to the upwind cell center,ρx,ρy,and ρzare the density gradient at the upwind cell center,which can also be computed by the K-exact least squares method mentioned before,qis the dynamic pressure,μ is the coefficients of the artificial viscosity defined as whereMais the Mach number in the upwind cell,Maca cutoff Mach number around 0.98,andCa constant normally between 1 and 2 used to increase the artificial viscosity.

3.3.Implicit algorithms and GMRES iteration

To achieve a rapid convergence of the solution,Eq.(5)is solved in an implicit scheme.Write Eq.(5)into an implicit fashion as

and the linearization of Eq.(13)results in a set of linear algebraic equations as

The linearized equations of Eq.(14)with a large sparseare solved by a generalized minimum residual(GMRES)iteration method with an incomplete lower upper(ILU)preconditioner.12The number of inner iterations in GMRES is chosen between 50 and 100 in the present study.The convergence criterion tol,which also significantly influences the solution efficiency,is suggested by a large number of numerical experiments as

With a converged iterative solution at each step,the potential of each cell is updated as

where λ is the relaxation factor used to control the rate of convergence and the stability.As the flow is relatively stable in a subsonic case,λ can be assigned as 1 from the beginning.However,in a transonic case,λ is suggested to start with a small positive value and then gradually increase to 1 when the residual decreases to a certain level.

3.4.Solution adaptation

Solution adaptation is applied to improve the accuracy of the flow solution.Due to the isentropic and irrotational assumptions of the full potential flow,velocity divergence is chosen as a criterion to refine/coarsen the grid cells.By doing so,the large gradient in the field with the flow changing dramatically,particularly near the shock wave,can be accurately captured.13Considering the influence of the cell size,the criterion parameter for refinement/coarsening is defined as where▽·Vis derived by the method of gradient reconstruction mentioned in Section 3.2.2,and the velocity on the interface is used for interpolation.liis the arithmetic cube root of the cell volume andr,an adjustable parameter,is assigned as 2 here.As a threshold,the root-mean-square of all τiis given as

wherenis the total number of grid cells.If τi＞ σd,cellishould be refined;if τi＜ 0.1σd,cellishould be coarsened.For a cell created by refinement,its potential φiis interpolated from the original cell as

where φjandVjare the potential and velocity of the original cell,respectively,andrjiis the vector from the centroid of the original cell to that of the new created cell.The potential of a cell created by coarsening is replaced by the average of the potentials of all coarsened cells.

4.Boundary conditions and viscous drag prediction

4.1.Wall boundary condition

Because of the inviscid and irrotational assumptions of the full potential equation,the wall boundary is non-penetrating.This type of boundary condition is implemented by a ghost cell method.12The potential value of a ghost cell can be given in a simple way as

However,the result computed by this simple method,as shown in Fig.5(a),indicates that this traditional method is not able to ensure the wall parallel velocity vector.To improve the accuracy of the wall boundary condition,this paper proposes an analytical method for assigning the value of φghost.

For the non-penetration condition,the velocity on the wall satisfies

The velocity on the wall can be written as the gradient of potential in the same way used in Section 3.2.2 for potential gradient reconstruction,that is,

whereCxis the coefficients of the relative φ in Eq.(10),so areCyandCz.By substituting Eq.(23)into Eq.(22),we have

Fig.5(b)represents the result computed by the analytical method.It is obvious that the velocity vector is parallel with the wall precisely,which indicates that a stricter wall boundary condition can be enhanced by the proposed method compared to the traditional simple way.

4.2.Kutta condition

Due to the irrotational assumption of the full potential flow,the Kutta condition should be properly applied so that a lift force can be produced.In the present study,the Kutta condition is implemented on cells and faces adjoining a smooth cut plane.Before generating a Cartesian mesh,this cut plane is firstly generated,starting from the trailing edge of an aircraft wing and extending to the far field along the trailing wake,and then is used to identify wake cells and wake faces during grid generation.It should be mentioned that the generating process of the cut plane and the subsequent Kutta condition implementation are fully automatic and thus do not attenuate the automation of the whole analysis.This procedure is summarized as follows:

(1)Find a vertex located on the trailing edge of a wing as a starting point.An approximate position near the trailing edge,e.g.,near the wingtip trailing edge,specif i ed beforehand by a user,is input along with an STL geometric model,as shown in Fig.6(a).The triangular facets close to this approximate position can be found via the aforementioned ADT search algorithm,and then a vertex on the trailing edge can be determined.

(2)Identify the whole trailing edge from the starting point by applying the ADT search algorithm again on the STL geometric model,as illustrated in Fig.6(a).

(3)Generate a cut plane from the trailing edge.According to the geometric information of the trailing edge,a cut plane is generated along the freestream direction.This generated plane is presented as an unstructured triangu-lar surface in an STL format,the same as the geometric model has,as illustrated in Fig.6(b).

(4)During generating and processing Cartesian grids,cells intersecting with the cut plane are divided so that the flow domain is separated by wake faces where the Kutta condition would be implemented,as illustrated in Fig.7.

After the above procedures,the circulation of each wake face should be calculated for Kutta condition implementation.As the circulation in each wing section varies along the wing span,a set number of points at the trailing edge are chosen to calculate the local circulation,which is given as

where Γ is the local circulation in each wing section,φtop,teand φbottom,teare the potential values on the top and bottom surfaces at the trailing edge,respectively.Both of them can be computed from the potential values of cells nearby by the K-exact least squares method mentioned in Section 3.2.The circulation at each wake face centroid is given by a linear interpolation from two circulations of selected adjacent points on the trailing edge,as shown in Fig.8.

Numerical experiments have shown that direct correction of the potential values of adjacent cells near the cut plane cannot ensure a convergence of the solution.The reason is that the cell size may vary greatly between these adjacent cells.To overcome this difficulty,a new method is proposed in this paper.To keep things simple,we take a 2D case as example.As illustrated in Fig.9,rather than directly correcting the potential values to implement the Kutta condition,the potential values of those cells in the stencils for gradient reconstruction of a wake faceABare corrected temporarily when the velocities of wake faces are calculated by

where φrightand φleftare the potential values of cells on both sides of the cut plane,φright,tempand φleft,tempare their temporary values,and ΓABis the circulation at the centroid of wake faceAB.By doing so,the potential gradient near the wake face is smoothed,and therefore,good convergence can be achieved.

4.3.Farfield boundary condition

The farfield boundary condition is implemented by correcting the potential values of cells adjacent to the farfield boundary.After computing the wing circulation in each step,the potential value of a farfield boundary cell is given as

whereV∞is the freestream velocity,Xthe coordinate of the cell center,θ the azimuthal angle measured from the Kutta face,and Γ is given by linear interpolation from wing circulation distribution by using the method shown in Fig.8.

4.4.Viscous drag prediction

The total drag consists of lift-induced drag,wave drag,and viscous drag.The induced drag and wave drag are solved by the full potential method in this paper,and the viscous drag is predicted by an engineering method.21The viscous drag coefficients for a single component is calculated by this engineering method,which is based on the boundary layer theory and the component shape factor correction method.The aircraft total viscous drag coefficients is estimated by a sum of all components’viscous drag coefficients.The geometric information for this engineering method could be obtained by the mesh module of the full potential solver in this paper.

5.Validations and results

5.1.Subsonic flow over DLR-F4

The first test case is a subsonic flow over a DLR-F4 wing-body configuration,and the conditions of the freestream areMa∞=0.6 and lift coefficientsCL=0.5.Computations were performed on Cartesian grids with solution adaptation,as shown in Fig.10.Note that the grids are refined in the field with the flow changing dramatically.Fig.11 shows the convergence histories,in which good convergence is achieved for all the solution adaptions.In the figure,the vertical axis means summing the absolute value of the residual results ofnsteps,and dividing the initial results,and taking the logarithm.As can be seen,the proposed Kutta condition method could provide satisfactory solution convergence.One thing should be mentioned that when the residual decreases to a certain amount,the circulation is frozen to reduce the remaining residual.On the same Cartesian grids,the computation time of this potential solver is decreased at least 75%compared with that of an Euler solver developed by our research group.That solver employs an implicit finite volume scheme with the GMRES algorithm to solve the Euler equations.The subsonic wing computation using this full potential approach on Cartesian grids consisting of 295412 cells requires about 4 min of wallclock time on a single processor of an Intel Xeon E5-2630 v2 CPU.

An illustration of the flow field is shown in Fig.12,which displays pressure coefficientsCpcontours.The surface is colored with pressure coefficients distributions.Fig.13 shows a comparison between the computed pressure coefficients distributions and experimental data22at seven semi-span sections,y/bmeans the relative spanwise position andx/cmeans the relative chordwise position.Generally speaking,the agreement between the numerical simulation and experimental data is fairly good.

5.2.Transonic flow over ONERA M6

The second example is a well-documented steady flow problem:transonic flow past an ONERA M6 wing.The flow conditions for this case areMa∞=0.839 and α =3.06°,and a Lambda-type shock wave is expected to appear on the upper surface of the wing.23Solution adaptation is applied to capture the shock wave accurately.As shown in Fig.14,grids are obviously refined at positions where the gradients of flow variables are large,especially at the shock wake position.Fig.15 shows the residual and lift coefficients convergence histories,and good convergence is achieved for all the solution adaptions.Without solution adaptation,the transonic computation on grids consisting of 151233 cells requires about 15 min of wall-clock time using the same processor in Section 5.1.

Fig.16 provides a quantitative assessment of the solution through pressure profiles at six spanwise positions.This figure displays the comparison of the pressure coefficients distribution between applying solution adaptation or not,which indicates that solution adaptation can improve accuracy.The inboard stations correctly display the double-shock on the upper surface,while these two shock waves merge together near a station at 0.8 semi-span to form a single strong shock wave in the outboard region of the wing.Fig.16 also gives the experimental data,23which show that although the solution locates this rear shock wave slightly behind its experimental data,a good agreement is observed.These discrepancies are attributed to the fact that this method is inviscid and ignores viscous effects in the boundary layer.Mach number contours on the wing surface and symmetrical plane are shown in Fig.17,which clearly shows a sharply captured Lambda-type 3D shock structure formed by the two inboard shock waves on the upper surface of the wing.

5.3.Computations on a BWB configuration

The purpose of the work in this paper is to analyze unconventional aircraft configurations.Therefore,for instance,an unmanned aerial vehicle(UAV)with a blended wing body(BWB)is used for simulation in subsonic and transonic flows.For the subsonic case,wind tunnel tests have been conducted for this BWB UAV in a low-speed wind tunnel at Nanjing University of Aeronautics and Astronautics,as shown in Fig.18.The initial computation grids before solution adaptation and after geometric adaptation are illustrated in Fig.7.

For the subsonic case, the flow conditions areMa∞=0.1763 and α =-6.28°to 9.42°.A grid convergence study is firstly considered for the simulation of α =-0.072°.As shown in Fig.19,the computed lift coefficients converges to a certain value as the grids are refined.The lift coefficients obtained by Richardson extrapolation based on the computed lift coefficients is 0.1847,and the difference between this value and the experimental data is within 2.5%.Fig.20 presents the comparisons between the computed and experimentalCL,CD,L/D,and polar.TakingCLas an example,the maximum error is less than 4%between this full potential method and the experimental data,which shows good accuracy of the method in this paper and indicates the capacity of this method in application for unconventional aircraft configuration studies.

For the transonic case,the flow conditions areMa∞=0.8 and α =0°.Fig.21 presents the Cartesian grids on the wall boundary after solution adaptation,in which the grids are refined by solution adaptation and the shock wave on the upper surface is caught successfully.Because of the lack of experimental data,we compare the results of the full potential solution with those obtained by the Euler solver developed by our research group on the same Cartesian grids.As shown in Fig.22,the surface pressure coefficients distributions demonstrate similar results of the transonic shock wave on the aftpart of the wing upper surface by the two methods,although the position of the shock wave computed by this full potential method is a little behind the Euler results,which is caused by the isentropic assumption.24The comparison indicates that the accuracy of the approach in this paper is satisfying.

6.Conclusion

A FVM for solving the full potential equation on adaptive Cartesian grids is proposed in this paper.The proposed method has been validated by simulating ONERA M6,DLR-F4,and a BWB UAV in subsonic or transonic flows.Conclusions are summarized as follows:

(1)Cartesian grids with geometric adaptation are firstly generated automatically with boundary cells processed by cell-cutting and cell-merging algorithms,and a common unstructured data structure is applied to store the topological information so that different kinds of complex CUT-cells could be described conveniently.

(2)The nonlinear full potential equation is discretized by a finite volume scheme on these Cartesian grids and iteratively solved in an implicit fashion with a GMRES algorithm.During computation,solution-based mesh adaptation is also applied so as to capture flow features more accurately.

(3)An improved ghost-cell method is proposed to implement the non-penetration wall boundary condition where the velocity-potential of a ghost cell is modified by an analytic method instead so that the wall boundary condition is better satisfied.

(4)According to the characteristics of the Cartesian grids,the Kutta condition is applied by specially computing the gradients on Kutta-faces without directly assigning the potential jump to cells adjacent wake faces,which can significantly improve the solution converging speed.

With the above-mentioned measures,the mesh can be generated automatically,and the computation time of this potential solver is decreased at least 75%compared with that of the Euler solver developed by our research group on the same Cartesian grids.Therefore,the proposed method can provide satisfying accuracy and rapid convergence in subsonic and transonic flow simulations with highly automatic grid treatment,indicating its potential application in aircraft conceptual design.

The next step in our development will be focusing upon viscid-inviscid coupling calculation and parallel computation.Future investigations will also study the application in highlift configuration simulation and other acceleration techniques.

Acknowledgements

This work was co-supported by the National Natural Science Foundation of China(No.11672133)and the Fundamental Research Funds for the Central Universities.The support from the Priority Academic Program Development(PAPD)of Jiangsu Higher Education Institutions is also appreciated.

1.Zhu ZQ,Wang XL,Wu ZC,Chen ZM.Multi-disciplinary optimization and numerical simulations in civil aircraft design.Acta Aeronaut et Astronaut Sin2007;28(1):1–13[Chinese].

2.de Weck O,Agte J,Sobieszczanski-Sobieski J,Arendsen P,Morris A,Spieck M.State-of-the-art and future trends in multidisciplinary design optimization.Reston(VA):AIAA;2007.Report No.:AIAA-2007-1905.

3.Holst TL.Transonic flow computations using nonlinear potential methods.Prog Aerosp Sci2000;36(1):1–61.

4.Jameson A.Iterative solution of transonic flows over airfoils and wings,including flows at Mach 1.Commun Pure Appl Math1974;27(3):283–309.

5.Johnson FT,Tinoco EN,Yu NJ.Thirty years of development and application of CFD at Boeing commercial airplanes,Seattle.Comput Fluids2005;34(10):1115–51.

6.Rubbert PE,Bussoletti JE,Johnson FT,Sidwell KW,Rowe WS,Samant SS,et al.A new approach to the solution of boundary value problems involving complex configurations.In:Noor A,editor.Computational mechanics—advancesamp;trends Amd.New York:ASME;1986.p.49–84.

7.Young DP,Melvin RG,Bieterman MB,Johnson FT,Samant SS,Bussoletti JE.A locally refined rectangular grid finite element method:application to computational fluid dynamics and computational physics.J Comput Phys1991;92(1):1–66.

8.Jameson A,Caughey DA.A finite volume method for transonic potential flow calculationsProceedings of AIAA 3rd computational fl uid dynamics conference.Reston(VA):AIAA;1997.p.35–54.

9.Caughey DA,Jameson A.Numerical calculation of transonic potential flow about wing-body combinations.AIAA J1979;17(2):175–81.

10.Caughey DA,Jameson A.Progress in finite-volume calculations for wing-fuselage combinations.AIAA J1980;18(11):1281–8.

11.Bolsunovsky AL,Buzoverya NP,Karas OV,Skomorokhov SI.An experience in aerodynamic design of transport aircraftProceedings of the 28th international congress of the aeronautical sciences.

12.Neel RE.Advances in computational fluid dynamics:turbulent separated flows and transonic potential flows[dissertation].Blacksburg(VA):Virginia Polytechnic Institute and State University;1997.

13.De Zeeuw DL.A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations[dissertation].Ann Arbor(MI):University of Michigan;1993.

14.Aftosmis MJ,Berger MJ,Melton JE.Robust and efficient Cartesian mesh generation for component-based geometry.AIAA J1998;36(6):952–60.

15.Bonet J,Peraire J.An alternating digital tree(ADT)algorithm for 3D geometric searching and intersection problems.Int J Numer Meth Eng1991;31(1):1–17.

16.Day MS,Colella P,Lijewski MJ,Rendleman CA,Marcus DL.Embedded boundary algorithms for solving the Poisson equation on complex domains.Berkeley(CA):Ernest Orlando Lawrence Berkeley National Laboratory;1998.Report No.:LBNL-41811.

17.Yang G,Causon DM,Ingram DM.Calculation of compressible flows about complex moving geometries using a three-dimensional Cartesian cut cell method.Int J Numer Meth Fluids2000;33(8):1121–51.

18.Sutherland IE,Hodgman GW.Reentrant polygon clipping.Commun ACM1974;17(1):32–42.

19.Aftosmis MJ,Berger MJ,Adomavicius G.A parallel multilevel method for adaptively refined Cartesian grids with embedded boundaries.Reston(VA):AIAA;2000.Report No.:AIAA-2000-0808.

20.Hafez M,South J,Murman E.Artificial compressibility methods for numerical solutions of transonic full potential equation.AIAA J1979;17(8):838–44.

21.Bonner E,Clever W,Dunn K.Aerodynamic preliminary analysis system II:Part I Theory.Washington,D.C.:NASA;1991.Report No.:CR-182076.

22.Redeker G.DLR-F4 wing-body configuration,a selection of experimental test cases for the validation of CFD codes.Paris:AGARD;1994.Report No.:AGARD AR-303.

23.Schmitt V,Charpin F.Pressure distributions on the ONERA-M6-wing at transonic Mach numbers,experimental data base for computer program assessment.Paris:Fluid Dynamics Panel Working Group 04;1979.Report No.:AGARD AR-138.

24.Mason WH.On the use of the potential flow model for aerodynamic design at transonic speeds.Reston(VA):AIAA;1995.Report No.:AIAA-1995-0741.

18 October 2016;revised 6 December 2016;accepted 22 December 2016