Project Topics
Get Complete Project Material File(s) Now! »
Jacobian computation and system resolution
The accuracy of the adjoint directly depends on the precision of the computed Jacobian and the inversion of the linear system. The Jacobian has thus to be computed as precisely as possible. In particular, it must not include the mass matrix. Additionally, in order to improve the di↵erentiation of inviscid fluxes, and contrary to the main solver where the first order Jacobian is considered, we perform MUSCL extrapolation on each edge before computing contributions to the Jacobian. Though, we do not di↵erentiate limiters as it may lead to too non-linear erratic e↵ects.
The absence of mass matrix and the use of a more precise di↵erentiation lead to a sti↵er system that SGS solver cannot handle. As we need the linear system to be converged by at least eight orders of magnitude (10−8) to obtain an accurate adjoint, we need a stronger linear solver. Thus we use a Krylov GMRES solver. It is still preconditioned with a LU-SGS solver, i.e. the vectors added to the Krylov subspace are not computed as AkJ but as the solution obtained with a single iteration of LU-SGS from the current residual.
Spalart-Allmaras Quadratic Constitutive Relation (SA-QCR)
Each turbulence model is known to fail at modelling some flows. In particular, the Boussinesq hypothesis, i.e. the assumed proportionality of the turbulent stress tensor and the strain tensor, is known to be false in some cases. For example, impinging jets or recirculating flows are classical counter example to this hypothesis. It is thus expected that turbulence models that relie on the Boussinesq hypothesis poorly predict such flows. In order to improve the prediction of recirculating flows that appear at the junction of the body and the wing of an aircraft, [Spalart 2000] proposed the following Quadratic Constitutive Relation.
The turbulent equations are the same as for the standard SA turbulence model and the Boussinesq hypothesis is amended so that the turbulent stress tensor depends on the shear stress and rotation tensor through a quadratic relation: TQCR = TSA − Ccr1 / O · T T SA + TSA · OT 0.
RANS solution initialization
The computation of RANS problems is significantly slower than inviscid flows. This is partly due to the boundary layer that is particularly slow to converge. When starting a computation from nothing, the choice of the initial field has thus a strong influence on the convergence of the problem. In order to improve this convergence, we use di↵erent initial fields. The idea is to compute at a very low expense an initial analytical field that will significantly improve the shape of the boundary layer. The solver needs at least an initial physical field that has a positive density and a positive pressure. The basic initialization consist in setting the whole domain to a reference constant flow. This obviously produce a very steep (inexistent) boundary layer that will require initially a low CFL and many iterations to propagate.
In order to ease this and be able to start with a higher CFL, we need to smooth the boundary layer. To do so we use a linear blending from the wall to the external flow. In order to further improve this initialization, we need to carefully chose the blending distance. Too small it will yield a too steep boundary layer. Too big it will produce the opposite problem, a boundary layer that takes a lot of time to get thinner.
To evaluate this optimal distance, we relie on boundary layer theory. The idea is that for blu↵ bodies, we will have a decent description of the growth of the boundary layer with a flat plate. This follows the same approach as for boundary layer mesh generation. We first identify a theoretical stagnation point, namely the leading point of the body. We can then use the empirical law giving the boundary layer thickness along the body /(x) ⇡ 0.16x (Rex)1/7 .
Periodic Adjoint
We also have to take into account the periodicity in the adjoint resolution. In the main solver we were able to simplify the matrix contribution by the use of ghost vertices. However, this rely on a direct modification of the linear system solver in Wolf. As we use another linear system solver in the adjoint resolution, using the same approach would require to apply the same modifications in the external library. Hence, we choose to properly take into account the periodicity in the matrix instead for adjoint resolution.
Figure 3.2 shows the periodic computational domain of a LS89. Vertex C is the periodic of A and ghost vertex D is the periodic of B. The periodicity hypothesis implies nonetheless that D has the same value as B but also that C has the same value as A. Hence, even though only D is generated as ghost vertex, C should be considered as well. Additionally, as the nodal values of D are copied from B, the residual of A actually depends on B.
Matrix construction: Contrary to the main solver, the full periodic matrix is assembled in the adjoint system. Periodic vertices are first identified (D is periodic of B, C is periodic of A, …) and the proper connectivity is reconstructed:
• Any periodic vertex is considered to be connected to no one: the line associated to C and D in the matrix are fixed to identity and the entries linking A to D and B to C are set tot 0 (not allowed).
• If an inner vertex is connected to a periodic vertex (like A and D), then it is reconnected to the proper inner periodic vertex (A is reconnected to B).
Non-linear Corrector computation
The numerical simulation of engineering problems involves a discrete solution which is alleged to converge toward the continuous solution of the problem as the size of the elements of the mesh decreases. As the exact analytical solution of the problem is sought, the di↵erence between the numerical and the analytical solution can be seen as a noise. In an engineering context, this noise can have disastrous consequences on the numerical prediction which can, in turn, lead to actual accidents or misconceptions.
It is thus essential to reduce and estimate this error. This is usually done by leading a grid convergence study, where the error is estimated with the di↵erence between two solutions computed on two successives grids. However, this approach is valid only if the asymptotic regime has already been reached and requires an additional finer grid (thus expensive) to appropriately estimate the error of the current solution. Moreover, the prescription and generation of a sequence of nested grids is time consuming and requires the expertise of an engineer. Mesh convergence study is mandatory to have a reasonable confidence into the discrete outputs, but it will eventually fail. Problems in mesh convergence are often related with insufficient resolution of local small details of the flow. Mesh adaptation techniques aim at adressing these two difficulties by automatically generating multiple successive grids following a sensor which estimates the error related to the local mesh size [Loseille 2007a, Alauzet 2010a]. It does not require any mesh prescription, but it needs a relatively accurate estimation of the numerical error to be efficient and to stop when the mesh and the solution are converged.
These estimations are based on the fact that discrete solution does not satisfy the continuous problem and the projection of the exact solution on the discrete space does not satisfy the discrete problem either. This can be stated in an a priori and an a posteriori ways : h(⇧hW) = « h 6= 0.
2D Turbulent NACA0012
The corrector is then validated on a turbulent case with the same geometry, at a Reynolds number of 5 ⇥ 105 with the Sapalart-Allmaras turbulence model. The Mach number is fixed to M = 0.7 and the angle of attack to ↵ = 1#. Adapted meshes are generated iteratively with feflo.a using a hessian based error estimate based on the Mach number field. An example of the resulting meshes and solutions is shown in Figure 4.23. Additionally, we generate for each mesh its h/2 subdivision and computed the solution on it for comparison in order to analyse the influence of mesh adaptation on the performance of the corrector.
Figure 4.25 show the convergence of the drag computed with a first and second order scheme for both the initial (red) and corrected (black) solution. The blue curve represents for each mesh the drag computed on its h/2 mesh subdivision which actual complexity is thus four times larger. Both curves exhibits the same trend as the previous laminar example, but it is now possible to quantify how mesh adaptation improves the solution on a mesh with four times more vertices compared to the actual h/2-mesh subdivison. This dims the apparent efficiency of the nonlinear corrector. In particular, we can see that the corrector is very efficient with a first order scheme as it yields the same solution as the h/2-mesh. It is slightly less efficient with a second order scheme but still satisfying.
Table of contents :
Introduction
I RANS solver implementation for turbulent applications.
1 Main solver description
1.1 Main flow solver description
1.1.1 Modeling equations
1.1.2 Spatial discretization
1.1.3 Boundary Conditions
1.1.4 Rotating Flow
1.1.5 Time integration
1.2 Linear System Resolution
1.2.1 Standard linear solvers
1.2.2 SGS – Line solver
1.2.3 Multigrid solver
1.2.4 Multigrid in Wolf:
1.2.5 Residual computation and limiter freezing
1.3 Adjoint
1.3.1 Functional sensitivity to the solution
1.3.2 Jacobian computation and system resolution
1.4 Verification and validation
1.4.1 Verification
1.4.2 Validation
2 RANS solver description
2.1 Turbulence modeling
2.1.1 Spalart-Allmaras Negative Model
2.1.2 Spalart-Allmaras Quadratic Constitutive Relation (SA-QCR)
2.2 Spatial discretization
2.2.1 Di↵erentiation
2.3 Wall Functions
2.3.1 RANS solution initialization
2.3.2 RANS adjoint
3 Periodic computations for turbomachinery applications
3.1 Implementation
3.2 Periodic Adjoint
3.3 Applications
3.3.1 LS89 case
3.3.2 RO37 case
4 Solution correction for error estimation
4.1 Non-linear Corrector computation
4.1.1 1D RANS Example
4.1.2 Formal discretization
4.1.3 Finite Volume Corrector
4.1.4 Finite Element Corrector
4.2 Numerical results: manufactured solutions
4.2.1 Finite elements 1D nonlinear advection di↵usion problem
4.2.2 Finite volume 2D Navier-Stokes equations
4.3 Numerical results: applied cases
4.3.1 2D Laminar NACA0012
4.3.2 2D Turbulent NACA0012
4.3.3 3D ONERA M6 wing turbulent Case
4.3.4 3D inviscid supersonic case: Sonic Boom Prediction Workshop .
4.4 Conclusion and perspectives
II Anisotropic mesh adaptation for RANS applications
5 Continous mesh framework for mesh adaptation
5.1 Continous mesh
5.1.1 A simple 1D example
5.2 Extension to higher dimensions
5.2.1 Euclidian and Riemannian metric spaces
5.2.2 Unit mesh
5.2.3 Operations on metrics
5.3 The continuous mesh framework
5.3.1 Duality discrete-continous: a new formalism
5.3.2 Continuous linear interpolation
5.3.3 Summary
5.4 Multiscale mesh adaptation
5.4.1 The non-linear adaptation loop
5.4.2 Optimal control of the interpolation error and optimal meshes
5.4.3 General error control
5.4.4 Error equidistribution
5.4.5 Application to numerical solutions
5.4.6 Control of the error in Lp norm
5.4.7 Conclusion
6 RANS mesh adaptation
6.1 RANS goal-oriented mesh adaptation
6.1.1 General non-linear error estimation
6.1.2 RANS goal-oriented error estimation: Linearization
6.1.3 RANS goal-oriented error estimation: Case of Spalart-Allmaras equations
6.1.4 Second derivatives error estimation
6.1.5 Hessian computation
6.2 Mesh adaptation cycle for applied cases
6.3 Applied cases: NASA Rotor 37
6.3.1 Geometry description
6.3.2 Mesh comparison
6.4 Applied cases: M6 wing
6.4.1 Case description
6.4.2 Feature-based mesh adaptation
6.4.3 Comparison of hessian reconstructions
6.4.4 Goal-oriented mesh adaptation
6.4.5 Comparison with other codes
6.5 Applied cases: Trap wing
6.5.1 Case description
6.5.2 Results
6.6 Applied cases: High-Lift Common Research Model
6.7 Improvement of RANS mesh adaptation
6.7.1 Adaptation with Wall-Laws: Subsonic NACA0012 validation
6.7.2 Hessian boundary correction: impact on a 2D High-Lift configuration
6.7.3 Conclusion
6.8 Conclusion
7 Periodic mesh adaptation
7.1 Periodic mesh adaptation
7.1.1 Handling the metric
7.1.2 Periodicity recovery
7.2 Application to turbomachinery
7.2.1 LS89 case
7.2.2 RO37 case
Conclusion
Conclusion and perspectives
Acknowledgments
Bibliography
