Adomian Decomposition Based Numerical Scheme for Flow Simulations

by   Imanol Garcia-Beristain, et al.

Time efficiency is one of the more critical concerns in computational fluid dynamics simulations of industrial applications. Extensive research has been conducted to improve the underlying numerical schemes to achieve time process reduction. Within this context, this paper presents a new time discretization method based on the Adomian decomposition technique for Euler equations. The obtained scheme is time-order adaptive; the order is automatically adjusted at each time step and over the space domain, leading to significant processing time reduction. The scheme is formulated in an appropriate recursive formula, and its efficiency is demonstrated through numerical tests by comparison to exact solutions and the popular Runge-Kutta-DG method.



page 1

page 2

page 3

page 4


Strang splitting schemes for N-level Bloch models

We extend to the N-level Bloch model the splitting scheme which use exac...

A Gauss-Seidel projection method with the minimal number of updates for stray field in micromagnetic simulations

Magnetization dynamics in magnetic materials is often modeled by the Lan...

Optimal-order error estimates for Euler discretization of high-index saddle dynamics

High-index saddle dynamics provides an effective means to compute the an...

ISALT: Inference-based schemes adaptive to large time-stepping for locally Lipschitz ergodic systems

Efficient simulation of SDEs is essential in many applications, particul...

Numerical Solution of Compressible Euler and Magnetohydrodynamic flow past an infinite cone

A numerical scheme is developed for systems of conservation laws on mani...

Doubly-Adaptive Artificial Compression Methods for Incompressible Flow

This report presents adaptive artificial compression methods in which th...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Computational fluid dynamics (CFD) is the numerical simulation of fluid flow, heat and mass transfer, turbulence, and other physical properties of fluids. This is achieved by solving Euler equations numerically for inviscid flows and Navier-Stokes equations for viscous flows. Details of the involved equations from modelling to numerical discretization can be found in (Versteeg and Malalasekera, 2007; Katate, 2013; Colonius and Lele, 2004b).

CFD is extensively used in a wide range of industrial sectors including aeronautics and automotive (Remaki L., 2011; Olander, 2011; Evans et al., 2011; Remaki et al., 2014), ventilation and gas turbine design (Califano and Steen, 2009; Moshfeghi et al., 2012; Axerio-Cilies and Iaccarino, 2012; Remaki et al., 2017; Wang1 et al., 2014; Yelmule and VSJ, 2013), chemical manufacturing (Singh et al., 2007), oil industry (Lu X and Xie P and Ingham D and Ma L, 2018), power generation and marine energy (Karthikeyan et al., 2016), food processing (Da-Wen, 2019; Anandharamakrishnan, 2013), water treatment (Samuelsberg and Hjertager, 1996), and others. An interesting overview of the application of CFD to industrial applications is provided in (Lutz, 2010; Ram et al., 2018). To illustrate this fact, one can mention, for instance, that wind tunnel testing in automotive and aeronautics sectors is strongly combined with CFD simulations for better and cost-effective results and bringing products to market faster. The exclusive use of CFD appears as a serious option in the near future. Bloodhound supersonic car fully designed with CFD is a good example of this trend(Remaki et al., 2014; Evans et al., 2010).

Aeroacoustics simulations, referred to by computational aeroacoustic (CAA), are among the most important CFD applications. The main objective of such simulations is noise reduction based on an accurate noise prediction. Reducing noise is of great importance for industry and human life quality. In transportation, a sensitive population’s mobility growth is expected thanks to fast transportation facilities’ development. In the US, the Joint Planning and Development Office (JPDO) is creating a new NextGen system plan for an expected air traffic increase by a factor of 3 towards 2025. Therefore reducing harmful sound effects becomes critical. A study of the impact of industrial noise on human life can be found in (Atmaca et al., 2005; Colonius and Lele, 2004b).

One of the most challenging aspects when dealing with a numerical discretization of Euler and Navier-Stokes equations is the time efficiency. Indeed in industrial applications, complex geometries are used for simulations resulting in big meshes size that significantly increases the processing time. The case of aeroacoustics applications is even worse; a very fine mesh, around 6 points per wavelength, is required, aggravated by the need for extended time propagation resulting in the use of large domains, adding to that the significant scale differences inherent to the nature of the phenomenon (Christopher, 1995).

CAA simulations can be classified into two categories:

direct noise computation (DNC) and hybrid methods (Colonius and Lele, 2004a). DNC method computes simultaneously acoustic and aerodynamic fields using compressible Navier-Stokes or Euler equations. The DNC requires low dispersive and low dissipative numerical methods, making the approach prohibitively expensive in terms of computational cost; consequently, its use limited for real-life applications (Colonius and Lele, 2004a). On the contrary, the hybrid method computes the aerodynamic and the acoustic fields in two separate steps, following Lighthill’s acoustic analogy (Lighthill, 1954, 1952). The approach requires first solving the aerodynamic field, used to evaluate an acoustical source that is then propagated. Linearised Euler equations (LEE) or nonlinear Euler equations (Liever et al., 2016), where the fluid is decomposed into a mean flow part and perturbation part, are the standard ways used to propagate the acoustic field. Most industrial CAA applications favor the hybrid method and often use LEE decomposition being faster and more robust than its nonlinear alternative (Liever et al., 2016; Schulze et al., 2017; Bissuel et al., 2018; Williamschen et al., 2015). Hybrid methods are significantly cost-effective compared to DNC, and are mostly adopted in aeronautics (Liever et al., 2016; Xiao-dong et al., 2015). Nevertheless, they still requiring high-order low-dispersion schemes, and hence long simulation time (Williamschen et al., 2015). In a conclusion, LEE is extensively used for aeroacoustic waves simulations due to the many advantages this approach offers (Schulze et al., 2017; Bissuel et al., 2018; Williamschen et al., 2015). Indeed, LEE in addition of being significantly faster than nonlinear Euler Equations, advection velocity, and amplitude, are easier to preserve by performing a linearization around the mean flow. In contrast, nonlinear Euler equations might completely damp waves by numerical diffusion. From another side, the information about the mean flow is preserved compared to a simple wave equation. Consequently, LEE offers a good compromise between simple wave propagation equations and nonlinear equations. Despite the linear nature of the LEE, due to the big required domains and very fine mesh, existing numerical methods need to be significantly improved especially for high-order methods.

Reducing computational time for high-order methods is a very active research field (Kroll et al., 2015). A significant effort is put into the implementation aspect using parallel architectures to attenuate the prohibitive time processing (Brown, 2010; Huerta et al., 2013). In parallel, considerable work is devoted to improving numerical schemes based on adaptivity techniques to avoid unnecessary computations.

Two main adaptive strategies are developed in literature; mesh adaptivity, referred to by -adaptivity, where the mesh quality is improved by performing local operations on mesh cells and based on a selected metric. The operations include refinement, swapping, and coarsening (DOLEJsi and FELCMAN, 2002; Peraire et al., 1992; Sørensen et al., 2003; Habashi et al., 2000; Remaki and Habashi, 2006). The other strategy, referred to by -adaptivity, consists of adapting the method order. Although order adaptivity is widely used in space (Zander et al., 2015; da Veiga et al., 2018), it has not been reported in time adaptivity, except adjustment of time-step size (Walter and Manera, 2016; Völcker et al., 2010). black

Within this context, this paper tries to open doors to build cost-effective time-adaptive (-adaptivity) numerical schemes through the semi-analytical Adomian decomposition technique for Euler Equations. Adomian decomposition method first proposed by Adomian Adomian and Rach (1983); Adomian (1994); Cherruault (1998), is used by many authors to solve a big range of problems covering linear and nonlinear PDEs either deterministic or stochastic. This demonstrates that the approach is a promising trend in the field of approximation of PDEs solutions. See, for instance, the important work of Wazwaz Ebaida et al. (2015); Duan et al. (2012); Wazwaz (2011); Rach et al. (2015); Singh and Wazwaz (2015); Momani and Odibat (2006); Khan et al. (2009); Momani and Odibat (2006), and other relevant articles Vahidi and Jalalvand (2012); Gbadamosi et al. (2012); Alabdullatif et al. (2007). The reader is particularly referred to the nice review paper by Duan et al Duan et al. (2012).

While almost all the work on Adomian decomposition in the literature focuses on finding semi-analytical solutions for PDEs (very limited work on CFD), we investigate the potential of using this technique to build efficient numerical schemes.

Very few works can be found on semi-analytical CFD solutions; we can cite solving Burgers’ equation (Mohyud-din et al., 2010; Zhu et al., 2010; Chen and An, 2008), and fractional derivative type problems such as the fractional Navier-Stokes equations (Khan et al., 2009; Momani and Odibat, 2006). This approach, however, is too prohibitive to be used for complex real-life problems. Indeed, the approach calculates successive Adomian series terms analytically, which become exponentially complex even for simple problems. Only the very few first terms can be actually calculated, resulting in a significant truncation error.

The use of Adomian decomposition to derive numerical schemes is even rarer; one can cite the work about the fractional Navier-Stokes equations (Birajdar, 2014). Here the method is applied to discrete Navier-Stokes. In other words, a Cartesian mesh is used for space discretization, which is a severe limitation for reel applications that uses complex CADs (computer-aided design) where unstructured meshes are required. Moreover the Riemann-Liouville integrals are used for time increment, which is numerically prohibitive. This makes the method not competitive compared to standard numerical methods.

The novelty of the present work can be summarized as follows; first, the proposed numerical method uses the same data structure used for standard numerical methods like finite element, finite volume, discontinuous Galerkin(DG) methods, and so on. Consequently, no extra complexity is required, and the method can be implemented in any existing CFD platform (commercial, open-source, or academic). Second, the method is time-exact, with no need for any integration formula. Finally, the method is time-adaptive which can sensibly improve the cost-effectiveness of existing space-discretization methods.

To derive the proposed method, the semi-analytical Adomian decomposition technique is applied to the Euler equations. By considering the decomposition at each time step, a time discretization scheme is obtained. A recursive formula is established and proved providing a simple and practical formulation of the scheme. We refer to the obtained method by ABS standing for Adomian Based Scheme. The most important property of the scheme is its local time-order adaptivity, in the sense that the order is automatically adjusted at each time step and over the whole space domain, conferring to the method its time-effectiveness.

Space discretization can be achieved by any standard method. In this work classical discontinuous Galerkin (DG) approach is selected, and ABS-DG is used to refer to the full discretization. Numerical tests are performed using ABS-DG and Runge-Kutta - DG (RK-DG), and then results are compared to demonstrate the performance of the proposed method. A comparison to exact solutions is also presented.

In section 2, a short overview of the Adomian decomposition method and the nonlinear Euler and linearized Euler equations (LEE) are given. In section 3, details of the proposed ABS and ABS-DG schemes are provided as well as a stability analysis. In section 4, the connection of ABS-DG to the Runge-Kutta (RK) in the linear case is established. Tests are performed and reported in section 5.

Ii Adomian Decomposition Overview

ii.i The Adomian Decomposition Method

In the following, a short description of the Adomian decomposition method is provided, for more details see Adomian (1994); Adomian and Rach (1983); Duan et al. (2012).

Given a differential equation, the method consists in decomposing the equation as follows,


where is the unknown, and are the linear parts of the differential operator, with being the easily invertible part. is the nonlinear part.

Note that it is not necessary to distinguish the non-invertible linear part from the nonlinear one. The sum can be represented by a single operator , reducing the decomposition to


then the solution and the operator are expanded as,




The Adomian polynomial coefficients are given by


and the terms of the series are computed as follows


ii.ii The non-conservative Euler equations

The two-dimensional compressible inviscid flow equations (Euler equations) formulated relative to a Cartesian coordinates system, in primitive variables, are given by


where and denote the averaged density and pressure of the fluid respectively and is the averaged velocity of the fluid in direction . Solutions to the set of equations are defined on a fixed spatial computational domain .

Euler equations can also be written in a vector form,


where, and


where is the sound speed.

ii.iii The Linearized Euler Equations (LEE)

LEE system is obtained by a linearization around a mean flow. This is achieved by decomposing the solution into mean and perturbation parts then substituting in equation (8) and assuming that


with and the matrices and evaluated at .

We obtain Rydin (2016); Imanol (2018); Blom (2003),


where and and are the matrices and evaluated at .

For a smooth background (mean flow) which implies that mean flow spatial derivatives are moderate, in the sense


equation (11) is simplified to


In this paper, matrices and are considered to be constant (a common simplification in aeroacoustics) without loss of generality. For convenience, the prime symbol is dropped from .

Iii Adomian Based Schemes (ABS)

In this section, a cost-effective numerical scheme to solve nonlinear Euler Equations is derived. The time discretization is obtained using the Adomian decomposition technique, while the DG method (any other discretization method can be used) is applied for the space discretization. We refer to the semi-discretized scheme by ABS standing for Adomian Based Scheme, and by ABS-DG to the fully-discretized scheme. The DG space discretization is implemented following the approach proposed by Shu Shu and Atkins (1998). For other classical DG discretizations, the reader is referred to Cockburns paper Cockburn et al. (1999).

iii.i ABS Scheme Derivation

To apply the Adomian decomposition technique described in section II.I to Euler equations (7), we propose to set

   and then    .

The nonlinear operator is given by the remaining terms,


with being the primitive variables


First is expanded as in (3)


the Adomian coefficients (4) are then give by


which yields

The series terms are computed recursively by (6),


with a four components vector


Let’s develop the first component , corresponding to the continuity equation.

From (17)

In the first summation part of , permute the order of the derivatives (assuming the required regularity to do so)

Using the Leibniz formula, we get

Note that

we then obtain


Similarly, the second summation part is simplified as,


Substituting (20), (21) into , we get the formula


A similar formula for is obtained by repeating the same steps.

For and , first expand as power series (note that can be considered close to zero since we are concerned by the limit)

that is,

A recursive formula is then obtained for

Substituted in (17) we get,

Using the same simplifications as for , we get the following formula for


In practice, time integration (18) is, in general, not easy to compute accurately. Indeed, Adomian coefficients are polynomials in time, we need to store all coefficients for an accurate time integration which can be a severe limitation for a large . The following Theorem remedy this problem, and an exact time integration is obtained by a simple multiplication by time .

Theorem 1: The Adomian coefficients can be expressed as


where , depends only on and .

Proof: We proceed by induction. The formula is proved for () and for the other components of the proof can be obtained in a similar way.

Note first that if property (24) is fulfilled by , it is fulfilled by as well. Indeed


Note also that the indices of the product terms in always sum to , this implies that satisfies the property (24) as long as satisfies it, that is


Now to initialize the induction process, let’s verify the property for and . From (23) the first Adomian coefficient


and then


The second Adomian coefficient is given by

and from (26)

Substituting and into , it clear that .

Now assume property (24) satisfied till index and let us prove it for . Since the relation is valid for index we have



Substitute the last expression in formula (23) we obtain

This shows that is a monomial of degree in time with a coefficient depending only on , that is

which achieves the proof of the Theorem. ∎

iii.ii The ABS formula

Using the previous theorem we can derive a formula for that doesn’t require any time integration, indeed


which leads using (23) to the following recursive formulation for that requires no time integration


Note that the formulation 31 provides a fully analytical scheme, depending only on the initial condition and its derivatives of all orders. To calculate analytically the ABS series terms (or Adomian series terms in general), the only way is to express each term in terms of leading to very complicated expressions even for the first few terms of the series. This is the case even for simple problems, as can be verified in Gbadamosi et al. (2012); Alabdullatif et al. (2007); Rach et al. (2015); Momani and Odibat (2006). Another handicap, the derivatives are defined in the strong sense, which requires the initial condition to be or at least with being the first truncated term in the Adomian series. Unfortunately, this is not suitable for real-life problems where solutions show in general discontinuities as in the CFD field, the target applications of this paper.

Alternatively, we propose to consider space derivatives in the weak sense and estimate them numerically. Although any numerical method can be used, in this work, the DG method is selected and implemented following the approach proposed by Shu

Shu and Atkins (1998). However, different DG method formulations can be found in the literature, see, for instance, Cockburn Cockburn et al. (1999). We refer to the obtained numerical scheme by ABS-DG.


The extension to is straightforward since the system of equations is the same with an additional equation for the third momentum component which is similar to other components. Extension to Navier-Stokes equations is straightforward as well, we can easily check the recursive propriety for viscous terms and get


iii.iii Time-adaptive property of the ABS

From eq. (25) we have

and the exact solution is given by the series,

The solution is a power series in time. In practice, the solution is approximated by the first first terms based on a given tolerance on the truncation error (), that is

with given by the series residual

shows that is the time-order.

It is clear that once the tolerance fixed, the first neglected term of the series (the order) will not be the same for each position and time step . Consequently, the order is automatically adjusted at each time step and domain position to the needed accuracy.

iii.iv The ABS-DG for LEE

For the LEE case, formula (31) wrote

iii.v Space discretization

Space discretization is achieved by applying the DG method for each term of the ABS series following Shu’s approach, see Shu and Atkins (1998) and (Lummer, 2016) for details. Blom’s extensive work on DG for LEE (Blom, 2003) is a good reference as well.

In the following, zero-order DG case, which corresponds to a finite volume scheme, is however presented, since it will be used subsequently for stability analysis.

Each term is approximated at the cell center by the average value

where is a given cell surface (for a two-dimensional domain). We then get

is the approximated solution, and corresponds to the index for which is smaller than a given tolerance.

Approximating fluxes as in the classical finite volume, the ABS-DG zero-order is obtained,


where fluxes can be approximated for instance using Lax-Friedrichs formulation


with indices referring to nodes and and to the values of the matrices at the corresponding nodes. Stability Analysis

The stability analysis is performed for the one-dimensional linear wave propagation equation,


For such a purpose, a zero-order ABS-DG in space is used. For convenience, the spatial discretization index is written as a superscript and ABS iteration as a subscript. This non-standard notation is intended to differentiate ABS iterations from classical finite difference time steps, where usually refers to the current time level. ABS-DG terms () using 32 are given by


Round-off errors () are also governed by the same equation.