3.  Work Programme

 Within this Project we aim at addressing a broad variety of problems that arise in the modeling of fluid dynamical processes.
In section 3.1 we report a list of such problems, while in section 3.2 we mention which kind of numerical methods will be considered for their numerical simulation.
Finally, in section 3.3 we provide details on several specific issues that will be addressed in this project. Although non exhaustive, the list of section 3.3 outlines some of the areas of major current interest in CFD; we discuss the type of investigation that will be carried out, as well as the kind of contribution that our project is expected to bring about for the solution of these problems.

 3.1 List of problems

 Our research activity will concern the mathematical modeling and analysis of fluid flows and it may include:

- Compressible and incompressible flows and design of models and algorithms able to cope with both cases at the same time;
- Single -and multi-phase fluid  flows;
- Reactive and turbulent fluid flows;
- Flows in porous media;
- Free-surface flows;
- Viscoelastic and non-Newtonian  fluid flows;
- Modeling of turbulent flows (two-equations models, Reynolds stress models, LES);
- Multifields mathematical and numerical models, in particular for the coupling of fluids and solids, viscous and inviscid flows, rotational and irrotational flows, molecular and continuous  regimes;

 Let us also mention a few examples of slightly more specialized problems that could be addressed:
- Electrochemistry and simplified models for electromagnetism;
- Advanced models for semiconductor devices;
- Shape design optimization in aero-dynamics.

We also wish to  point out that we insist on such an impressive list of important problems  both from a scientific and an industrial viewpoint, since progress on one of these topics will certainly lead to some progress on the others.

3.2 List of methods

 Both classical and modern numerical  methods , such as:
- finite differences, finite volumes and finite elements,
- spectral and high order methods,
- wavelet and  multi resolution analysis,
- multiscale methods, including multigrid and domain decomposition  methods,
 will be addressed.

 Assessment of grid quality and control, especially for complex geometries, analysis of a-priori and a-posteriori error estimates, development of parallel preconditioners, coupling of heterogeneous numerical techniques, represent additional issues of interest for our research program.

3.3 Specific fields to be addressed

 We provide hereby some additional details for a few instances of  fields that we aim at investigating in our project.

3.3.1 Computation of low speed flows on structured grids

 Because of greater ease in grid generation, a main trend in computational fluid dynamics is currently the development of computing methods using unstructured grids. However, for very computing time and storage intensive application, structured grids are more attractive. Examples are direct and large-eddy simulation of turbulence. This has been feasible until now only on very simple grids (orthogonal and quasi-uniform), and is mostly done on Cartesian grids. Extension to more general geometries and grids can in the near future only be foreseen for structured boundary-fitted grids. For incompressible flows (Mach < 0.15, say) or flows containing low-speed regions an additional difficulty crops up. A choice has to be made between staggered and nonstaggered grids. Current methods for compressible flow break down for small Mach number, although this can be alleviated somewhat by preconditioning (Briley, Choi, Merkle, van Leer, Turkel, Viviand ....) or introduction of artificial compressibility (Guerra, Gustafsson, Merkle, Metcalfe ...). These measures usually falsify transient behaviour, thus ruling out unsteady applications. But they use nonstaggered grids, making discretization on general structured and unstructured grids relatively easy. Truly incompressible methods obviously do not break down at low Mach numbers, but the classical dependable staggered grid discretization is harder to implement accurately on general structured grids than nonstaggered discretizations; even more so on unstructured grids. Nevertheless, progress has been made in this direction (Davidson and Hedberg, Ikohagi nad Shin, Kwak, Rosenfeld and Vinokur, Shyy, Wesseling, ...). Time dependence poses no special problem in this approach. Nonstaggered schemes need to be used in conjunction with special measures, such as pressure weighted interpolation (Rhie and Chow, Demirdzic and Peric, Leschziner, Rodi, ...). Here troubles arise with time dependence, as for the artificial compressibility approach (Chorin, ...).
 This whole area needs to be developed further in order to transfer direct and large-eddy simulation of turbulence to more general geometries. Furthermore, such further development is also necessary to bridge the gap between the incompressible and the compressible domain in a better way than is possible at present. This is required for technological applications such as low speed flows with large density variations (for example in combustion or study of catalysts in chemistry) or flows incorporating both low speed and high speed zones (flows in inlets of internal combustion engines, or in chemical reactors with high-speed gas injectors, studies of pressure vessel rupturing, and even, as has been recently realized, accurate treatment of compressible wing flows near the leading edge stagnation line).
 To this must be added improvement of algorithmic efficiency of domain decomposition for greater geometric flexibility, and parallel computing methods, using both parallelism over subdomains and inside subdomains. The mathematical background of domain decomposition has seen much progress over the past ten years, but, unavoidably, there is a gap between the class of equations for which theory has been developed and the equations of fluid dynamics.

3.3.2 Multidimensional Upwinding Techniques in CFD

 Developed initially for the approximation of steady state solutions of the two-dimensional Euler equations on unstructured grids, multidimensional upwind methods comprise a decomposition step, in which an appropriate discrete form of the system of equations is split into simple components, and a distribution step, in which each component is allowed to evolve in time via a technique in which a cell-based quantity, known as the fluctuation, is distributed to the nodes of the grid.  These genuinely multidimensional upwind schemes are based on the concept of fluctuation distribution where the underlying representation of the solution is a continuous piecewise linear function with the unknowns stored at the grid nodes.  This differs radically from other finite volume algorithms and in some ways more closely resemble a finite element method, although it is generally more complicated to introduce the concept of upwinding into finite element schemes.  The need for truly multidimensional schemes for the solution of the scalar advection equation has long been recognised by many researchers as intrinsic to the solution of multidimensional hyperbolic systems, and not only in conjunction with multidimensional upwinding.  In the case of the multidimensional flux balance distribution methods appropriate scalar schemes began to appear about ten years ago and a number of upwind schemes now exist which are both positive (monotone in one dimension) and linearity preserving (roughly second order accurate at the steady state).  These are now well understood, particularly after recent work which has drawn a close parallel between the higher order two-dimensional flux limited schemes.  The focus of most current research on these schemes is the decomposition stage.  The decomposed forms of the Euler equations can roughly be divided into two groups, simple wave models and approximate diagonalisations.  The decomposition of the system into components corresponding to simple wave solutions of the equations was first proposed by Roe and has led to a wide variety of suggestions such as those presented in the related steady pattern model.  The characteristic decomposition of Deconinck et al which appeared at roughly the same time as the first simple wave models, has recently been modified, first by modifying the transformation used in the partial diagonalisation of the system and then by attempting to diagonalise a preconditioned form of the equation.  This final group of preconditioned methods has at last provided the basis for schemes which can be used reliably to obtain accurate steady state solutions to the two-dimensional Euler equations for a variety of configurations.

3.3.3 Mesh Adaptation in Computational Applied Mathematics

 The quality of the solutions obtained using multidimensional upwinding methods, or any other numerical scheme for the solution of partial differential equations, depends on the grid upon which the discretisation is defined.
Generally, it can be expected that the solution quality will be improved in a particular region of the flow if the nodes are more highly concentrated there.  Moreover, regions of the flow where the solution varies rapidly are likely to require a higher density of nodes than areas of relatively uniform flow to achieve a similar degree of accuracy in the approximation.  Therefore, if the grid is restricted to have a certain number of nodes, the solution quality depends on the distribution of the nodes within the domain.
 This indicates that if the equations are to be solved efficiently as well as accurately then it is worthwhile attempting to distribute the nodes of the grid in an intelligent manner, one which will in some way depend on the solution itself.  In other words the initial grid usually produced with the aid of only intuitive knowledge of the solution, should be adapted in some way to the numerical approximation to the solution.
 In two dimensions, the usefulness of grid adaptation can be illustrated simply by considering the modeling of air flow around an isolated airfoil section.  An aerodynamicist requires the fast and accurate prediction of the lift and drag created by the aerofoil, two quantities which are affected strongly by the position and strength of any shocks in the flow, where the gradient of the solution is very high. In a shock capturing method such as multidimensional upwinding nothing is known about the position of the shocks a priori and so a good grid (one which the solution error is low relative to the number of nodes in the grid) can only be obtained by applying an adaptation algorithm once information about the solution is known, so that nodes are concentrated around the shocks.
 Although the optimal grid is unattainable due to lack of an exact error measure or, for that matter, exact information about the solution, it may be approached by attempting to distribute some local error measure evenly over the grid. This can be done using an algorithm which moves grid points towards regions where the error in the solution is likely to be high, e.g. where the flow gradient is high.
It is more common to attempt to reduce the error by selectively refining the mesh to increase the node concentration in regions of the flow where the chosen error measure is high. Also, on unstructured triangular grids, the error may be reduced by  edge swapping, since many quasi one-dimensional flow features, such as shocks, are captured more accurately when they are aligned with grid.
However, both of these alternatives are to some degree achieved by careful movement of the nodes, a technique which has the advantage of changing neither the number of nodes in the grid or their connectivity. In fact, if the node movement algorithm is inexpensive and the pseudo time stepping algorithm with which it is combined is robust, the only additional expense arises from the usual decrease in the allowable time-step produced by smaller grid cells, a problem that is partially alleviated by the use of local time-stepping.  The ultimate goal is to produce a grid on which the global error of the discrete solution is minimised for the available number of nodes.

3.3.4 Multifields (heterogeneous) modeling in computational fluid dynamics

 The numerical simulation of complex flows of real life interest is very often a hard task. The computational efficiency of the adopted numerical method can be enhanced by resorting to different strategies. Among these, the so-called multifield approach can be regarded as a way for simplifying both the mathematical description of the problem at hand and the associated numerical algorithm, via the combination of multiple kinds of differential models and numerical schemes. The resulting approach has a "heterogeneous" nature that can allow a substantial reduction of the computational efficiency without affecting the accuracy of the computed solution.
 Our project will address some interesting instances of multifield methods for the modeling of flow problems. In particular we will consider the case of external flows in aerodynamics. When high Reynolds numbers are considered, it is known that apart from a thin region around the body surface, the velocities of fluid flows around streamlined bodies are of the same order of magnitude as the free-stream velocity. In our approach we will proceed to the coupling between Navier-Stokes and Euler equations through a fictitious interface which can be drawn just beyond the boundary layer (Glowinski, Periaux, Quarteroni, Canuto, Brezzi, Russo, Pironneau,..).
 The mathematical analysis of this multifield model enlightens the correct matching conditions that the physical variables should satisfy across the interface: continuity of fluxes, and continuity of those Riemann invariants associated to characteristic lines pointing into the far field (Euler) region. Based on a suitable splitting of these interface conditions, an effective method can be set up to solve the Euler problem independently of the Navier-Stokes one, and conversely.
 This approach can also be adopted for flows where transition between laminar and turbulent regimes occur at fixed places. In such cases the flow-field can be split in three regions: the outer region where Euler equations are used, a laminar flow region upstream the transition points described by the Navier-Stokes equations, and a turbulent flow region described by the Reynolds-averaged Navier-Stokes equations together with  additional transport equations for the turbulent kinetic energy and its dissipation rate.
 Another example that we aim at is the coupling of instationary compressible Navier-Stokes equations modeling the flow around airfoils or three-dimensional wings, with the compressible stationery hyperbolic Euler equations further away which can be modeled by the Poincaré decomposition of the velocity and a combination of the vorticity transport equation with a divergence-free velocity and the potential flow satisfying a potential equation with the density governed by a pressure function. The far-field finally is modelled by a stationary flow including the vorticity transport. The aim is the correct modeling of the transmission conditions between the different regimes and the formulation and exploitation of nonlocal boundary conditions taking into account the far-filed behaviour. (Engquist, Pironneau, Quarteroni, Wendland, ...)
 The numerical methods that are used for approximating the various models can be chosen according to the nature of the subproblem at hand. In particular, shock capturing methods and centered high order methods can be alternated in the proper subregions, through a transferring mechanism at the interface which can also account for a  geometrical nonconformity.
 The time advancement scheme can also be handled in a separate manner, for instance using local time steppings obeying different CFL conditions on the different flow regions, especially when approximating steady flows.
 Similar situations are encountered in the simulation of  the full potential equation for transonic flow calculation (Wendland, Coclici, Hafez,...) where, for instance, boundary elements can be coupled with finite elements.
 Other instances of multifield numerical modeling are provided by fractional step methods that treat the different terms of the overall differential problem sequentially, giving rise to independent subproblems that can be faced by numerical methods of different nature. An example is provided by the solution of the Navier-Stokes equations for incompressible flows by the so-called projection method. Another instance arises from the simulation of multiphase flows in porous media (e.g. in oil reservoirs, Espedal, Ewing, Chavent, Jaffré, Quarteroni,...).
 Another topic that we aim at addressing is the use of low-order methods as preconditioners for high order methods, in connection with the solution of flow problems by domain decomposition strategies (Orszag, Deville, Mund,...). This approach combines the inherent accuracy enjoyed by polynomial methods with the geometric flexibility and low complexity of finite element methods.
At this stage, domain decomposition algorithms can provide valuable computational tools within parallel computing environments.

3.3.5 Stabilization of numerical methods in CFD

The application of the Galerkin finite element method to certain problems of mathematical physics and engineering (including CFD) yields numerical approximations that are deficient in that they do not inherit the stability properties of the continuous problem.
 Stabilized methods constitute a systematic methodology for improving stability behavior without compromising accuracy. As such, stabilized methods have provided fundamental solutions to the problem of discrete approximations in several practically important areas.
 One of the difficulties involved with stabilized methods like SUPG is the dependence of the schemes by certain parameters, whose value is not completely determined by the theory. The numerical solution is quite sensitive to the fine tuning of these parameters, so that the robustness of the algorithm can be poor.
 In recent years several efforts have been undertaken to better understand the theoretical foundations and origins of stabilized methods.
 One of the most promising approach is based on the identification of stabilized methods with Galerkin finite element methods employing finite element spaces enriched with so-called "bubble" functions.
 This observation may be found in several works, (Pierre, Brezzi, Bristeau, Franca, Mallet, Roge, Arnold, et al). Subsequently it was realized that strict equivalence of stabilized methods and bubble function methods could only be accomplished with specially constructed bubbles (Baiocchi, Brezzi, Franca). The most recent rendition of this philosophy derives the bubbles from certain element-level boundary-value problems. These are referred to as "residual-free bubbles" (Brezzi, Russo, Franca). They have provided a satisfactory basis for the derivation of stabilized methods with excellent properties and employ well established and accepted concepts throughout their construction. The conceptual viewpoint is to attack the original problem first with Galerkin's method involving standard and simple polynomial finite element spaces (e.g., linear triangles) and then improve the approximation, and correct any deficiencies with regard to stability, by systematically enriching the space with residual-free bubbles. Alternatively, the method can be viewed as simply the Galerkin finite element method applied to the enlarged finite element space consisting of the standard polynomials and residual-free bubbles. This observation obviates any possible criticism of the resulting stabilized method because its construction emanates from such a universally accepted derivational procedure, namely, Galerkin's method for a particular finite-dimensional subspace enjoying good properties for the problem under consideration.
 The resulting scheme is completely parameter-free, in the sense that for instance in the case when the SUPG method is reproduced, the stabilization constants can be explicitly computed.
 The residual-free bubble component of the solution can be also used as an "a posteriori" error indicator working in all regimes (convection-dominated and diffusion-dominated) (Russo).
 Another approach which has also provided a firm theoretical foundation for the construction of stabilized methods is the "variational multiscale method" presented recently (Hughes). The philosophy of this method is different from the previous one. The view here is that the standard Galerkin finite element method with simple polynomial spaces is an inadequate numerical paradigm for many practically important problems, in particular, those involving fine scale features that are numerically unresolvable due to the length scale of elements composing the mesh. It is acknowledged in this view that even if one is uninterested in resolving, or "seeing", the fine scale features, their effect on the coarse, or resolvable, scales must be accounted for in order to accurately calculate the coarse scales. Thus the implementation of the variational multiscale procedure consists of two steps: The first is purely non-numerical - the original problem is decomposed into two subproblems. One involves solving for the fine scales in terms of the coarse scales. The result is substituted into the second subproblem which results in a modified problem involving only the coarse scales. This is sometimes referred to as a "subgrid-scale model". This problem turns out to be a suitable one for presentation to the standard Galerkin finite element method employing simple polynomial-based elements. Because unresolved scales have been removed, a successful approximation naturally follows. Application of the Galerkin method to the modified problem is the second and final step of the variational multiscale approach. Various practical approximations can be made within the variational multiscale procedure. For example, one can assume that the fine scale phenomena exists only in the interior of element subdomains, before introduction of, and regardless of the particular approximating finite element spaces to be introduced in the Galerkin step. This assumption leads to a framework which permits an identification with stabilized methods and, of particular interest, an equivalence with residual-free bubbles (Brezzi, Franca, Hughes, Russo).
 The identification of residual-free bubbles and variational multiscale method is an important step toward a clarification of the stabilization procedures for the finite element method.
 In the framework of the present project, we plan to apply the methodologies described above to finite element approximations of flow problems. This work should help to develop robust and parameter-free stabilization procedures to be readily used in the applications.