Circular Function-Based Gas-Kinetic Scheme for Simulation of Viscous Compressible Flows

Abstract

A stable gas-kinetic scheme based on circular function is proposed for simulation of viscous compressible flows in this paper. The main idea of this scheme is to simplify the integral domain of Maxwellian distribution function over the phase velocity and phase energy to modified Maxwellian function, which will integrate over the phase velocity only. Then the modified Maxwellian function can be degenerated to a circular function with the assumption that all particles are distributed on a circle. Firstly, the RAE2822 airfoil is simulated to validate the accuracy of this scheme. Then the nose part of an aerospace plane model is studied to prove the potential of this scheme in industrial application. Simulation results show that the method presented in this paper has a good computational accuracy and stability.

1 Introduction

With the development of numerical simulation, the computational fluid dynamics (CFD) is becoming more and more important in industrial design of aircraft since its high-fidelity description of the flow field compared to the engineering method. Most numerical schemes are based on directly solving the Euler or N-S equations [1,2,3]. Whereas, a new method we proposed here is to solve the continuous Boltzmann model at the micro level and N-S equations at macro level. The gas-kinetic scheme (GKS) is commonly used as a continuous Boltzmann model which is based on the solution of Boltzmann equation and Maxwellian distribution function [4,5,6,7]. The GKS attracts more researchers’ attention during the last thirty years since its good accuracy and efficiency in solving the inviscid and viscous fluxes respectively [8, 9].

The GKS is developed from the equilibrium flux method (EFM) [10] to solve inviscid flows in the very beginning. Then, the Kinetic Flux Vector Splitting (KFVS) [11] scheme is applied to solve collisionless Boltzmann equation. The Bhatnagar Gross Krook (BGK) gas kinetic scheme was developed by Prendergast et al. [12], Chae et al. [13], Xu [14, 15] and other researchers based on the KFVS scheme. The particle collisions are considered in BGK scheme to improve the accuracy, which contributes great developments and application potential for BGK gas kinetic scheme.

In this work, a stable gas-kinetic scheme based on circular function framework is proposed for simulating the 2-D viscous compressible flows. Most existing GKS are based on Maxwellian function and make it time consuming and complex. Hence, the original Maxwellian function, which is the function of phase velocity and phase energy, is simplified into a function including phase velocity only. The effect of phase energy is contained in the particle inter energy $e_{p}$. Then, based on the assumption which all particles are concentrated on a circle, the simplified Maxwellian function can be reduced to a circular function, which makes the original infinite integral to be integrated along the circle. Compressible flow around RAE2822 airfoil is simulated to validate the proposed scheme. Furthermore, the nose part of an aerospace plane model is studied to prove the application potential of this scheme.

2 Methodology

2.1 Maxwellian Distribution Function

Maxwellian distribution function is an equilibrium state distribution of Boltzmann function. The continuous Boltzmann equation based on Bhatnagar Gross Krook (BGK) without external force collision model is shown as Eq. (1):

$$ \frac{\partial f}{\partial t} + \xi \cdot \nabla f = \frac{1}{\tau }\left( {f^{eq} - f} \right), $$

(1)

Where $f$ is the gas distribution function and the superscript eq means the equilibrium state approached by $f$ through particle collisions within a collision time scale $\tau$.

The Maxwellian distribution function is

$$ f^{eq} = g_{M} = \rho (\frac{\lambda }{\pi })^{{\frac{D + K}{2}}} e^{{ - \lambda \left[ {\sum\limits_{i = 1}^{D} {\left( {\xi_{i} - U_{i} } \right)^{2} + \sum\limits_{j = 1}^{K} {\zeta_{j}^{2} } } } \right]}} , $$

(2)

in which $U_{i}$ is the macroscopic flow velocity in i-direction and $\lambda = m/(2kT) = 1/(2RT)$. The number of phase energy variables is $K{ = }3 - D + N$. $D$ is the dimension and N is freedom rotational degree number.

The heat ratio $\gamma$ can be expressed as:

$$ \gamma = \frac{b + 2}{b} = \frac{K + D + 2}{K + D}, $$

(3)

in which $b$ represents the freedom degree number of molecules.

Based on Maxwellian function (Eq. (2)), the continuous Boltzmann equation (Eq. (1)) can be recovered to N-S equations by applying Chapman-Enskog expansion analysis [4] with following conservation moments equations:

$$ \int {g_{M} } d\Xi = \rho , $$

(4a)

$$ \int {g_{M} \xi_{\alpha } } d\Xi = \rho u_{\alpha } , $$

(4b)

$$ \int {g_{M} (\xi_{\alpha } \xi_{\alpha } + \sum\limits_{j = 1}^{K} {\zeta_{j}^{2} } )d\Xi = \rho (u_{\alpha } u_{\alpha } + bRT)} , $$

(4c)

$$ \int {g_{M} \xi_{\alpha } \xi_{\beta } d\Xi = \rho u_{\alpha } u_{\beta } + p\delta_{\alpha \beta } } , $$

(4d)

$$ \int {g_{M} (\xi_{\alpha } \xi_{\alpha } + \sum\limits_{j = 1}^{K} {\zeta_{j}^{2} } )\xi_{\beta } d\Xi = \rho [u_{\alpha } u_{\alpha } + (b + 2)RT]u_{\beta } ,} $$

(4e)

in which $d\Xi = d\xi_{\alpha } d\xi_{\beta } d\xi_{\chi } d\zeta_{1} d\zeta_{2} \cdots d\zeta_{K}$ is the volume element in the phase velocity and energy space. $\rho$ is the density of mean flow, the integral domain for each variable is $( - \infty , + \infty )$.

Due to phase velocity is independent from phase energy space, Eq. (2) can be written as:

$$ g_{M} = g_{M1} \cdot g_{M2} , $$

(5)

$$ g_{M1} = \rho \left( {\frac{\lambda }{\pi }} \right)^{{\frac{D}{2}}} e^{{^{{ - \lambda \sum\limits_{i = 1}^{D} {(\xi_{i} - U_{i} )^{2} } }} }} , $$

(6)

$$ g_{M2} = \rho \left( {\frac{\lambda }{\pi }} \right)^{{\frac{K}{2}}} e^{{^{{ - \lambda \sum\limits_{j = 1}^{K} {\zeta_{j}^{2} } }} }} . $$

(7)

If we define $d\Xi _{1} = d\xi_{\alpha } d\xi_{\beta } d\xi_{\chi }$ and $d\Xi _{2} = d\zeta_{1} d\zeta_{2} \cdots d\zeta_{K}$, then we can get $d\Xi = d\Xi _{1} d\Xi _{2}$. With these definitions, the integral form of $g_{M2}$ can be concluded as:

$$ \int {g_{M2} } d\Xi _{2} = \int {\left( {\frac{\lambda }{\pi }} \right)^{{\frac{K}{2}}} e^{{^{{ - \lambda \sum\limits_{j = 1}^{K} {\xi_{j}^{2} } }} }} d\zeta_{1} d\zeta_{2} \cdots d\zeta_{K} = 1.} $$

(8)

Substituting Eqs. (5)–(8) to Eq. (4), we have

$$ \int {g_{M1} } d\Xi _{1} = \rho , $$

(9a)

$$ \int {g_{M1} \xi_{\alpha } } d\Xi _{1} = \rho u_{\alpha } , $$

(9b)

$$ \int {g_{M1} (\xi_{\alpha } \xi_{\alpha } + 2e_{p} )d\Xi _{1} = \rho (u_{\alpha } u_{\alpha } + bRT)} , $$

(9c)

$$ \int {g_{M1} \xi_{\alpha } \xi_{\beta } d\Xi _{1} = \rho u_{\alpha } u_{\beta } + p\delta_{\alpha \beta } } , $$

(9d)

$$ \int {g_{M1} (\xi_{\alpha } \xi_{\alpha } + 2e_{p} )\xi_{\beta } d\Xi _{1} = \rho [u_{\alpha } u_{\alpha } + (b + 2)RT]u_{\beta } } , $$

(9e)

in which $e_{p}$ is particle potential energy, shown as Eq. (10):

$$ e_{p} = \frac{1}{2}\int {g_{M2} \sum\limits_{j = 1}^{K} {\zeta_{j}^{2} } } d\Xi _{2} = \frac{K}{4\lambda } = [1 - \frac{D}{2}(\gamma - 1)]e, $$

(10)

where $e = p/[(\gamma - 1)\rho ]$ is the potential energy of mean flow. It can be seen from Eqs. (9) and (10) that $e_{p}$ is independent from phase velocity $\xi_{i}$.

2.2 Simplified Circular Function

Suppose that all the particles in the phase velocity space are concentrated on a circle which has center $(u_{1} ,\,u_{2} )$ and radius $c$, shown as:

$$ \left( {\xi_{1} - u_{1} } \right)^{2} { + }\left( {\xi_{2} - u_{2} } \right)^{2} { = }c^{2} , $$

(11)

in which $c^{2}$ means the mean particle kinetic energy and we have

$$ c^{2} = \frac{{\int {\sum\limits_{i = 1}^{D} {(\xi_{i} - U_{i} )^{2} g_{M1} d\Xi _{1} } } }}{{\int {g_{M1} d\Xi _{1} } }} = \frac{D}{2\lambda } = D(\gamma - 1)e. $$

(12)

By substituting Eq. (12) into Eq. (9a) we can get Eq. (13), which is the mass conservation form in the cylindrical coordinate system:

$$ \rho = \int_{0}^{2\pi } {\int_{0}^{\infty } {\rho \left( {\frac{\lambda }{\pi }} \right)^{{\frac{D}{2}}} e^{{^{{ - \lambda r^{2} }} }} \delta (r - c)rdrd\theta } } = 2\pi c\rho \left( {\frac{\lambda }{\pi }} \right)^{{\frac{D}{2}}} e^{{^{{ - \lambda c^{2} }} }} . $$

(13)

Then we can get the simplified circular function shown as follows:

$$ g_{C} = \left\{ {\begin{array}{*{20}l} {\frac{\rho }{2\pi }} \hfill & {(\xi_{1} - u_{1} )^{2} + (\xi_{2} - u_{2} )^{2} = c^{2} ,} \hfill \\ 0 \hfill & {else.} \hfill \\ \end{array} } \right.. $$

(14)

All particles are concentrated on the circle and the velocity distribution is shown as Fig. 1. Then the phase velocity components in the Cartesian coordinate system can be expressed as:

Fig. 1.

Configuration of the phase velocity at a cell interface

$$ \xi_{1} = u_{1} + c\,\text{cos}(\theta ), $$

(15a)

$$ \xi_{2} = u_{2} + c\,\text{sin}(\theta ). $$

(15b)

Substituting Eqs. (14) and (15) to Eq. (9), the conservation forms of moments to recover N-S equations can be expressed as follows:

$$ \int {g_{M1} } d\Xi _{1} = \int_{0}^{2\pi } {g_{C} } d\theta = \rho , $$

(16a)

$$ \int {g_{M1} \xi_{\alpha } } d\Xi _{1} = \int_{0}^{2\pi } {g_{C} } \xi_{\alpha } d\theta = \rho u_{\alpha } , $$

(16b)

$$ \int {g_{M1} (\xi_{\alpha } \xi_{\alpha } + 2e_{p} )d\Xi _{1} } = \int_{0}^{2\pi } {g_{C} (\xi_{\alpha } \xi_{\alpha } + 2e_{p} )d\theta } = \rho (u_{\alpha } u_{\alpha } + bRT), $$

(16c)

$$ \int {g_{M1} \xi_{\alpha } \xi_{\beta } } d\Xi _{1} = \int_{0}^{2\pi } {g_{C} \xi_{\alpha } \xi_{\beta } } d\theta = \rho u_{\alpha } u_{\beta } + p\delta_{\alpha \beta } , $$

(16d)

$$ \int {g_{M1} (\xi_{\alpha } \xi_{\alpha } + 2e_{p} )\xi_{\beta } d\Xi _{1} } = \int_{0}^{2\pi } {g_{C} (\xi_{\alpha } \xi_{\alpha } + 2e_{p} )\xi_{\beta } } d\theta = \rho [u_{\alpha } u_{\alpha } + (b + 2)RT]u_{\beta } . $$

(16e)

2.3 Governing Equations Discretized by Finite Volume Method

N-S equation discretized by finite volume method in 2-dimensional can be expressed as:

$$ \frac{{d{\mathbf{W}}_{I} }}{dt} = - \frac{1}{{\Omega _{I} }}\sum\limits_{j = 1}^{{N_{f} }} {F_{Nj} } S_{j} , $$

(17)

in which $I$ represents the index of a control volume, $\Omega _{I}$ is the volume, $N_{f}$ means the number of interfaces of control volume $I$ and $S_{j}$ is the area of the interface in this volume. The conservative variables ${\mathbf{W}}$ and convective flux ${\mathbf{F}}_{n}$ can be expressed as:

$$ {\mathbf{W}} = \left[ {\begin{array}{*{20}c} \rho \\ {\rho u_{1} } \\ {\rho u_{2} } \\ {\rho (u_{1}^{2} + u_{2}^{2} + bRT)/2} \\ \end{array} } \right],\quad {\mathbf{F}} = \left[ {\begin{array}{*{20}c} {\rho U_{n} } \\ {\rho u_{1} U_{n} + n_{x} p} \\ {\rho u_{2} U_{n} + n_{y} p} \\ {U_{n}^{{}} \rho (u_{1}^{2} + u_{2}^{2} + (b + 2)RT)/2} \\ \end{array} } \right]. $$

(18)

Suppose that cell interface is located on $\varvec{r} = 0$, then the distribution function at cell interface (Eq. (19)) consists two parts, the equilibrium part $f_{{}}^{eq}$ and the non-equilibrium part $f_{{}}^{neq}$:

$$ f(0,t) = g_{C} (0,t) + f^{neq} (0,t) = f^{eq} (0,t) + f^{neq} (0,t). $$

(19)

To recover N-S equations by Boltzmann equation from Chapman-Enskog analysis [4, 9, 16,17,18,19], the non-equilibrium part $f_{i}^{neq} (0,t)$ applying Taylor series expansion in time and physics space can be written as:

$$ f^{neq} (0,t) = - \tau_{0} \left[ {g_{C} (0,t) - g_{C} ( - \xi \delta t,t - \delta t) + o(\delta t^{2} )} \right]. $$

(20)

Substituting Eq. (20) into Eq. (19) and omit the high order error item, then we have

$$ f(0,t) \approx g_{C} (0,t) + \tau_{0} \left[ {g_{C} ( - \xi \delta t,t - \delta t) - g_{C} (0,t)} \right], $$

(21)

in which $g_{C} ( - \xi \delta t,t - \delta t)$ is the distribution function around the cell interface, $\tau_{0} = \tau /\delta t = \mu /p\delta t$ is dimensionless collision time. $\mu$ is dynamic coefficient of viscosity, $\delta t$ is the streaming time step which represents the physical viscous of N-S equations. Now convective flux at cell interface can be expressed as:

$$ {\mathbf{F}} = {\mathbf{F}}^{\rm I} + \tau_{0} ({\mathbf{F}}^{{{\rm I}{\rm I}}} {\mathbf{ - F}}^{\rm I} ), $$

(22)

in which ${\mathbf{F}}^{\text{I}}$ represents contribution of equilibrium distribution function $g_{C} (0,t)$ at cell interface and ${\mathbf{F}}^{\text{II}}$ means equilibrium distribution function $g_{C} ( - \xi \delta t,t - \delta t)$ at surrounding point of the cell interface. Their functions are shown as:

$$ {\mathbf{F}}^{\text{I}} = \left[ {\begin{array}{*{20}c} {\rho u_{1} } \\ {\rho u_{1} u_{1} + p} \\ {\rho u_{1} u_{2} } \\ {(\rho E + p)u_{1} } \\ \end{array} } \right]^{face} , $$

(23)

$$ {\mathbf{F}}^{{{\rm I}{\rm I}}} = \int {\xi_{1}^{cir} } \varphi_{a}^{cir} g_{C}^{cir} d\theta , $$

(24)

in which the superscript face represents value at $(0,t)$ and cir means $( - \xi \delta t,t - \delta t)$.

When $\xi_{1}^{{}} \ge 0$, we have expressions as Eq. (25). The superscript L become R when $\xi_{1}^{{}} < 0$:

$$ \xi_{1}^{cir} = u_{1}^{L} - \frac{{\partial u_{1}^{L} }}{{\partial x_{1} }}(u_{1}^{ + } + c^{ + } \,\text{cos}(\theta ))\delta t - \frac{{\partial u_{1}^{L} }}{{\partial x_{2} }}(u_{2}^{ + } + c^{ + } \,\text{sin}(\theta ))\delta t + c^{ + } \,\text{cos}(\theta ), $$

(25a)

$$ \xi_{2}^{cir} = u_{2}^{L} - \frac{{\partial u_{2}^{L} }}{{\partial x_{1} }}(u_{1}^{ + } + c^{ + } \,\text{cos}(\theta ))\delta t - \frac{{\partial u_{2}^{L} }}{{\partial x_{2} }}(u_{2}^{ + } + c^{ + } \,\text{sin}(\theta ))\delta t + c^{ + } \,\text{sin}(\theta ), $$

(25b)

$$ e_{p}^{cir} = e_{p}^{L} - \frac{{\partial e_{p}^{L} }}{{\partial x_{1} }}(u_{1}^{ + } + c^{ + } \,\text{cos}(\theta ))\delta t - \frac{{\partial e_{p}^{L} }}{{\partial x_{2} }}(u_{2}^{ + } + c^{ + } \,\text{sin}(\theta ))\delta t, $$

(25c)

$$ g_{C}^{cir} = g_{C}^{L} - \frac{{\partial g_{C}^{L} }}{{\partial x_{1} }}(u_{1}^{ + } + c^{ + } \,\text{cos}(\theta ))\delta t - \frac{{\partial g_{C}^{L} }}{{\partial x_{2} }}(u_{2}^{ + } + c^{ + } \,\text{sin}(\theta ))\delta t. $$

(25d)

Substituting Eq. (25) into Eq. (24), the final expression of the equilibrium distribution function at the surrounding points of the cell interface ${\mathbf{F}}^{\text{II}}$ can be calculated by

$$ {\mathbf{F}}^{\text{II}} (1) = \int {\xi_{1}^{cir} } g_{C}^{cir} d\theta = \int {(a_{0} + a_{1} r + a_{2} s)} (g_{0} + g_{1} r + g_{2} s)d\theta , $$

(26a)

$$ {\mathbf{F}}^{\text{II}} (2) = \int {\xi_{1}^{cir} \xi_{1}^{cir} } g_{C}^{cir} d\theta = \int {(a_{0} + a_{1} r + a_{2} s)} (a_{0} + a_{1} r + a_{2} s)(g_{0} + g_{1} r + g_{2} s)d\theta , $$

(26b)

$$ {\mathbf{F}}^{\text{II}} (3) = \int {\xi_{1}^{cir} \xi_{2}^{cir} } g_{C}^{cir} d\theta = \int {(a_{0} + a_{1} r + a_{2} s)(b_{0} + b_{1} r + b_{2} s)(g_{0} + g_{1} r + g_{2} s)d\theta ,} $$

(26c)

$$ \begin{aligned} {\mathbf{F}}^{\text{II}} (4) = & \int {\xi_{1}^{cir} (\frac{1}{2}\left| {\varvec{\xi}^{cir} } \right|^{2} + e_{p}^{cir} )} g_{C}^{cir} d\theta \\ = & \int {(a_{0} + a_{1} r + a_{2} s)\left( {\frac{1}{2}\left( {(a_{0} + a_{1} r + a_{2} s)^{2} + (b_{0} + b_{1} r + b_{2} s)^{2} } \right) + (e_{0} + e_{1} r + e_{2} s)} \right)} \\ & (g_{0} + g_{1} r + g_{2} s)d\theta . \\ \end{aligned} $$

(26d)

The operators in Eq. (26) can be expressed as:

$$ \begin{aligned} a_{0} = u_{1}^{L} - \frac{{\partial u_{1}^{L} }}{{\partial x_{1} }}u_{1}^{ + } \delta t - \frac{{\partial u_{1}^{L} }}{{\partial x_{2} }}u_{2}^{ + } \delta t,a_{1} = c^{ + } - \frac{{\partial u_{1}^{L} }}{{\partial x_{1} }}c^{ + } \delta t,a_{2} = - \frac{{\partial u_{1}^{L} }}{{\partial x_{2} }}c^{ + } \delta t, \\ b_{0} = u_{2}^{L} - \frac{{\partial u_{2}^{L} }}{{\partial x_{1} }}u_{1}^{ + } \delta t - \frac{{\partial u_{2}^{L} }}{{\partial x_{2} }}u_{2}^{ + } \delta t,b_{1} = - \frac{{\partial u_{2}^{L} }}{{\partial x_{1} }}c^{ + } \delta t,b_{2} = c^{ + } - \frac{{\partial u_{2}^{L} }}{{\partial x_{1} }}c^{ + } \delta t, \\ e_{0} = e_{p}^{L} - \frac{{\partial e_{p}^{L} }}{{\partial x_{1} }}u_{1}^{ + } \delta t - \frac{{\partial e_{p}^{L} }}{{\partial x_{2} }}u_{2}^{ + } \delta t,e_{1} = - \frac{{\partial e_{p}^{L} }}{{\partial x_{1} }}c^{ + } \delta t,e_{2} = - \frac{{\partial e_{p}^{L} }}{{\partial x_{2} }}c^{ + } \delta t, \\ g_{0} = g_{C}^{L} - \frac{{\partial g_{C}^{L} }}{{\partial x_{1} }}u_{1}^{ + } \delta t - \frac{{\partial g_{C}^{L} }}{{\partial x_{2} }}u_{2}^{ + } \delta t,g_{1} = - \frac{{\partial g_{C}^{L} }}{{\partial x_{1} }}c^{ + } \delta t,e_{2} = - \frac{{\partial g_{C}^{L} }}{{\partial x_{2} }}c^{ + } \delta t, \\ \end{aligned} $$

(27)

in which $u_{1}^{ + } ,\,u_{2}^{ + }$ and $c_{{}}^{ + }$ are predictional normal velocity, tangential velocity and particle specific velocity at cell interface. These velocities can be obtained both by Roe average [19] and the value of the former moment at cell interface.

3 Numerical Simulations

To validate the proposed circular function-based gas-kinetic scheme for simulation of viscous compressible flows, the RAE2822 airfoil and nose part of aerospace plane model are discussed.

3.1 Case1: Compressible Flow Around RAE2822 Airfoil

The transonic flow around RAE2822 airfoil is discussed to validate the accuracy of the proposed scheme. The free stream has Mach number of $Ma_{\infty } = 0.729$ with 2.31° angle of attack.

Figure 2 shows the pressure contours obtained by the proposed scheme. Figure 3 shows the comparison between experimental data [20], Roe scheme and the presented scheme. It can be seen clearly that the presented scheme has good accordance with both experimental data and results of Roe scheme.

Fig. 2.

Pressure contours of presented scheme around RAE2822 airfoil.

Fig. 3.

Pressure coefficient comparison of RAE2822 airfoil surface.

3.2 Case2: Supersonic Flow Around Nose Part of Aerospace Plane Model

The supersonic flow around the nose part of an aerospace plane model (wind tunnel) is simulated in this section. The model is 290 mm in length, 58 mm in width, the head radius is 15 mm and semi-cone angle is 20°. The free stream has Mach number of $Ma_{\infty } = 3.6$ with −5, 0 and 5° angle of attack respectively. Figure 4 gives out the pressure contours simulated by the presented scheme at 5° angle of attack, and Fig. 5 is the Mach number contours at the same condition.

Fig. 4.

Pressure contours at 5° angle of attack

Fig. 5.

Mach number contours at 5° angle of attack

Figure 6 shows comparison of pressure coefficient distribution on upper and lower surface of aerospace plane model at −5 angle of attack, which computed by scheme presented, Roe scheme and Van Leer scheme. It can be seen that the result of present solver has good accordance with other two numerical schemes. Similar conclusions can be obtained according to Figs. 7 and 8. Therefore, the solver presented in this article shows both high computational accuracy and numerical stability. Hence, the circular function-based gas-kinetic scheme shows the potential of future industrial application in flight vehicle research and design.

Fig. 6.

Pressure coefficient comparison at −5° angle of attack

Fig. 7.

Pressure coefficient comparison at 0° angle of attack

Fig. 8.

Pressure coefficient comparison at 5° angle of attack

4 Conclusions

This paper presents a stable gas-kinetic scheme based on circular function. It is focus on improving the calculation efficiency of existing GKS. Firstly, simplifying the original Maxwellian function, which is the function of phase velocity and phase energy, into the function of phase velocity. Furthermore, reducing the simplified function to a circular function. Hence, the original infinite integral can be changed to the integral along the circle.

Transonic flow around RAE2822 airfoil and supersonic flow around the nose part of an aerospace plane model are studied. The results show a good computational accuracy and the potential for future industrial application.