A GAS KINETIC SCHEME FOR HYBRID SIMULATION OF PARTIALLY RAREFIED FLOWS

Approaches to predict §ow ¦elds that display rarefaction e ̈ects incur a cost in computational time and memory considerably higher than methods commonly employed for continuum §ows. For this reason, to simulate §ow ¦elds where continuum and rare¦ed regimes coexist, hybrid techniques have been introduced. In the present work, analytically de¦ned gas-kinetic schemes based on the Shakhov and Rykov models for monoatomic and diatomic gas §ows, respectively, are proposed and evaluated with the aim to be used in the context of hybrid simulations. This should reduce the region where more expensive methods are needed by extending the validity of the continuum formulation. Moreover, since for high-speed rare¦ed gas §ows it is necessary to take into account the nonequilibrium among the internal degrees of freedom, the extension of the approach to employ diatomic gas models including rotational relaxation process is a mandatory ¦rst step towards realistic simulations. Compared to previous works of Xu and coworkers, the presented scheme is de¦ned directly on the basis of kinetic models which involve a Prandtl number correction. Moreover, the methods are de¦ned fully analytically instead of making use of Taylor expansion for the evaluation of the required derivatives. The scheme has been tested for various test cases and Mach numbers proving to produce reliable predictions in agreement with other approaches for near-continuum §ows. Finally, the performance of the scheme, in terms of memory and computational time, compared to discrete velocity methods makes it a compelling alternative in place of more complex methods for hybrid simulations of weakly rare¦ed §ows.

p ρT (nondimensional equilibrium pressure) p r ρT r (nondimensional rotational pressure) p t ρT t (nondimensional translational pressure) q x Total heat §ux in x-direction q r x Rotational heat §ux in x-direction q t x translational heat §ux in x-direction R gas constant

INTRODUCTION
At intermediate altitudes (7090 km), the §ow around hypersonic aircraft can be characterized as mainly continuum with localized areas (generated by the rapid expansion in the wake of the vehicle as well as by strong gradients in shock waves and boundary layers) that display rarefaction e¨ects.The §ow conditions near the vehicle surface and in the wake determine the drag and the heat transferred to the vehicle and its payload.Therefore, it is important that these regions are simulated using appropriate physical models.When the gradients of the macroscopic variables become so steep that their length scale is of the same order as the average distance travelled by molecules between collisions, the number of impacts is not enough to drive the §uid towards a local thermodynamic equilibrium.At these conditions, the §ow can no longer be considered a continuum and the transport terms in the NavierStokes (NS) equations fail since the constitutive relation is not valid.
The mathematical model at molecular level is the Boltzmann transport equation (BTE) [1] and for the regions of the §ow ¦eld where highly nonequilibrium e¨ects occur, the direct simulation Monte Carlo (DSMC) method [2] is typically employed to statistically estimate the solution of the BTE.Alternatively, a discrete velocity method (DVM) [3,4] can be used to solve a kinetic model approximation of the BTE [57].
In previous works [811], hybrid techniques have been introduced to simulate §ow ¦elds where continuum and rare¦ed regimes coexist.In these methods, the more expensive approaches, such as DSMC or DVM, are employed only where needed and is coupled with a ¦nite-volume scheme for the NS equations used where the §ow is continuum.A hybrid technique couples two simulation methods by means of information exchange between the parts of the §ow domain.In recent works, this has been achieved using an overlap region where §ow state variables or numerical §uxes are exchanged between the two models [12] or employing a bu¨er region where the two models are blended at equation level [13].Recent works on these methods focused on rare¦ed high-Mach §ow can be found in [1416].
An alternative approach is the uni¦ed gas-kinetic scheme (UGKS) [17,18] which uses a ¦nite-volume method where the numerical §uxes are based on the solution of the Shakhov model [6] for a monoatomic gas or the Rykov model [7] for a diatomic gas with rotational nonequilibrium.Where the §ow is underresolved, by accounting for the pressure jump in the de¦nition of the collision time, additional numerical viscosity is added resulting in a shock thickness of the order of the cell size [19].This allows the UGKS to simulate §ows in both rare¦ed and continuum regimes.
In the present work, two GKS methods, analytically-de¦ned on the basis of the ChapmanEnskog (CE) expansion of nondimensional Shakhov and Rykov models, are proposed for the simulation of weakly rare¦ed §ows.The derivatives of the equilibrium function and the time derivatives of the primitive variables are analytically de¦ned, employing the compatibility condition of the kinetic model for the latter.In previous works from Xu and coworkers [20,21], similar GKS are de¦ned using the CE solution of the BhatnaganGrossKrook (BGK) model with rotational nonequilibrium and a scaling of the energy numerical §ux [19] to correct the Prandtl number.Moreover, in those schemes, the required derivatives are expressed in terms of Taylor series where the coe©cient are calculated by means of properties of the employed BGK model.The proposed GKS, due to the use of the CE expansion, is limited to near-continuum regions but is simpler than the UGKS [17,18].However, the validity of the approach can be extended considering a modi¦ed collision time [22].Also, this correction, in the present work, is de¦ned fully analytically for both schemes.
Based on a literature survey of related works, the authors believe that the proposed GKS represents an e©cient method, relative to DVM, capable of modeling complex diatomic gas §ows with moderate rarefaction e¨ects but with signi¦cant rotational nonequilibrium.As such, the proposed approach is a novel alternative to the DVM for a range of practically relevant §ows.Moreover, the update of the nonequilibrium distribution function, as used in the UGKS, is neglected reducing the memory cost of the approach.Thus, the GKS method represents a viable option to reduce the cost of hybrid simulations by reducing the region where the expensive method is needed and extending the validity of the continuum formulation.
The schemes are built in a computational framework described in section 2 that also includes a DVM for the kinetic Boltzmann equations successfully employed for di¨erent monoatomic cases [23].The framework has been recently improved [24] with the addition of the Rykov model and an Ellipsoidal-Statistical (ES) model [25] for diatomic gases with rotational nonequilibrium.In section 3, the mathematical de¦nition of the Rykov model and of the proposed GKS are described.Finally, in section 4, some preliminary results where the presented GKS is coupled with a DVM are shown and an assessment of reliability and performance of the GKS is presented.

MULTIPHYSICS CODE
The methods used in the present work are built in the multiphysics code (MC) developed at the University of Liverpool [23,26,27].Multiphysics code is a com-putational framework designed for simulations of complex §ows, where di¨erent mathematical §ow models are employed for di¨erent regions of the §ow domain depending on the §ow physics.The NS equations represent the baseline level of mathematical models used.For the continuum §ow solver based on the compressible NS equations, a cell-centered block-structured ¦nite-volume method is employed using the AUSM + /up (advection upstream splitting method) for the convective §uxes [28].For low-speed §ow analysis, the framework further includes a Lattice Boltzmann Method as well as a Vortex-In-Cell method for vortex-dominated incompressible §ows.
In the present work, emphasis is placed on the simulation of hypersonic, partially-rare¦ed §ows for which mathematical models at a more detailed level of physics than the compressible NS equations are required.For §ows with strong nonequilibrium and rare¦ed e¨ects, the framework includes Molecular Dynamics (MD) methods as well as deterministic DVM for a range of kinetic Boltzmann equations.The Shakhov and ES models are included for monoatomic gas §ows, while the Rykov model and a polyatomic ES-BGK model were implemented for diatomic gas §ow simulations [24].The kinetic models are discretized using a discrete velocity method within a ¦nite-volume method for multiblock structured grids; second order total variation diminishing (TVD) time marching is employed.The velocity (phase) space is discretized using either a uniformly spaced method with the trapezoidal rule for the evaluation of the moments of the distribution functions or a Gauss-quadrature method with modi¦ed Hermite polynomials.For the present high-speed §ow cases, the uniform velocity space with the trapezoidal rule is the preferred approach and was used exclusively.
The memory and CPU time requirements are considerable and for this reason, an e©cient parallel implementation involving ¢two£ levels of parallelism was conceived.In this parallelization, the phase space as well as the §ow domain are distributed over the processes.First, the phase space is partitioned in regular subspaces, each to be assigned to separate processes within separate message passing interface (MPI) communicators.The overall number of processes is then divided by the number of partitions to obtain the required number of communicators.The mesh-blocks in physical space are then distributed over these communicators to obtain an equal distribution of the cells.An important factor in the performance is the number of processes assigned to the velocity space discretization.Limiting the size of these communicators will minimize the overhead in collective operations required for the evaluations of the moments of the distribution functions.However, this implies an increasing number of communicators over which the mesh-blocks are divided, potentially creating a load imbalance in physical space.The best performance is obtained if the number of velocity space partitions is chosen such that the load imbalance in physical space is typically less than 10%.
The parallel performance of the coupled continuum/kinetic solver was investigated as part of a European Union funded project (PRACE Preparatory  4), µ is determined by Eq. ( 6), and specular wall boundary condition Access), providing access to the SuperMUC computer at the Leibniz Rechenzentrum (LRZ) in Munich.In the project, 2048 to 8192 cores of Intel R Xeon R processors were used for a range of test cases.The example shown here is for the Mach 8 §ow of a diatomic gas around the waverider geometry (leading-edge diameter 5λ and body length 5 • 10 3 λ) shown in Fig. 1a.A multiblock mesh with 867 blocks and 3 • 10 6 cells was used.The multiphysics simulations involved a region around the leading-edges which were simulated at the kinetic level; in the case shown in Fig. 1b, the kinetic- §ow domain comprises 300 blocks with 575,000 cells.A reduced spatial extent with about half the blocks was also considered, corresponding to the ¢thin£ layer results in Fig. 2. Velocity space discretizations as 24 3 and 32 3 were used divided in 8 × 4 × 4 or 8 × 8 × 4 partitions.The §ow demonstrates signi¦cant thermodynamic nonequilibrium in the kinetic- §ow domain around the leading-edges as can be seen from the rotational temperature relaxation in Fig. 1c.The strong scaling for these cases are shown in Fig. 2, showing a 85 percent parallel e©ciency from 2048 to 8192 processes.
For the coupling between the continuum solver and the kinetic solver, MC employs state-based or §ux-based hybrid techniques.In the state-based coupling (Fig. 3a), over an overlap region, the macroscopic variables are obtained from the microscopic solution while the latter is constructed from the macroscopic state.In the §ux-based coupling, instead, at a designed interface, the numerical §uxes for one solver are obtained from the other as shown in Fig. 3b, while for hybrid NS-DSMC simulations, the choice of a state-based coupling is the preferable one due to the lower scattering error that it involves [12].In the current deterministic DVM kinetic solver, the statistical scatter is absent, creating more §exibility in the used coupling technique.In the literature, the domain decomposition is generally done during the simulation on the basis of a continuum breakdown parameter, for example, as done by other researchers in [12,29].At the moment, this feature is still under development in MC and the DVM domain has to be de¦ned by the user in an input ¦le.However, the framework can perform a recon¦guration of the di¨erent regions throughout the calculation when the domain de¦nition ¦le is modi¦ed by the user.

Nondimensional Rykov Model
Considering the §ow of a diatomic gas, let us assume that the gas temperature is not too high, so that the vibrational degrees of freedom are not excited, and not too low, so that the rotational degrees of freedom can be considered to be fully excited.In this case, the particle distribution function f (x, c, t, ζ), which describes the state of the gas, will be a function not only of the spatial coordinate x, the particle velocity c, and the time t, but also of the rotational degrees of freedom ζ.The Rykov model represents an extension of the Shakhov model where also rotational nonequilibrium is considered and has been proved to be a reliable kinetic approximation, up to the heat §uxes moments of the BTE, for this kind of §ow [7,18,30,31].Since the rotational degrees of freedom are considered fully excited, ζ is reduced by the model and the second distribution function is obtained.
Employing the following nondimensional variables: the nondimensional distribution functions of the model result thus, for the Rykov model written in terms of F = mf , one obtains ; where and the total collision time τ is expressed as µ t /p t with the viscosity determined from the translational temperature.
In a system of colliding particles, energy is transferred between the various internal modes.These collisions tend to drive the internal energy distributions towards their respective equilibrium state and the number of them necessary to push a particular mode to the equilibrium is the collision number, Z, associated to that mode [32].The Rykov model is based on the assumption that the fraction of collisions involving the excitation of the rotational degrees of freedom, Z r , is a given constant or a function of the §ow temperatures.Several works provide an expression of Z r as a function of the temperature in the §ow ¦eld.Probably, the ¦rst attempt to appear in the literature is the theoretical work of Parker [33] where, employing an empirical nonimpulsive model and assuming coplanar collisions and zero initial rotational energy, the following approximate expression is obtained: where T * = 91.5 K is the characteristic temperature of the intermolecular potential and (Z r ) ∞ = 23.5 is the limiting value suggested in [34].While Parker£s expression ( 2) is derived involving a large number of simplifying assumptions, the overall dependence on the temperature is in agreement with the more rigorous treatment of [35].However, this expression does not involve any dependence on the di¨erent translational and rotational temperatures.Thus, in the recent literature, models derived from data ¦tting, either from numerical or experimental results, have been employed.In [7,30,31], the following expression for the collision number is presented to be used with the Rykov model: where ψ( " T ) = 0.767 + 0.233 " An alternative expression for Z r (T t , T r ) derived from MD simulations can be found in [36]: where a 1 = 1.33868; a 2 = −6.19992; a 3 = −0.00107942;and 0 < b ≤ 1.It is important to notice that considering the moments of the Rykov model collision term, the relaxation process in the model is described as ρ(T − T r )/(Z r τ ) while in [33,36], Jeans equation is considered leading to ρ(T t − T r )/(Z r τ ).This means that the collision number in the Rykov model results in For the viscosity law, Rykov and his coworkers [7,30,31] suggest otherwise, a simpler power law with an exponential factor of 0.72 [37] can be employed.To make the system (1) complete, the values of the constants δ, ω 0 , and ω 1 need to be determined.In [38], ω 0 = 0.2354 and ω 1 = 0.3049 or ω 0 = 0.5 and ω 1 = 0.286 are given for diatomic gases.Both pairs of values have been successfully employed in [30,31,38,39] with δ −1 = 1.55.In the present work, the values ω 0 = 0.5 and ω 1 = 0.286 are employed.
The dimensionless macroscopic quantities can be obtained from F 0 and F 1 by means of the following formulae:

Gas-Kinetic Scheme for Near-Continuum Flows Based on the Rykov Model
Integrating in time the nondimensional reduced Rykov model system (1) and taking the moments , the following can be obtained for the update of the nondimensional macroscopic variables: where the source term is and the velocity (phase) space is discretized using a uniformly spaced method with the trapezoidal rule for the evaluation of the moments of the distribution functions.As previously shown in [40,41], reconstructing the time-dependent distribution functions at the cell-faces, i. e., F 0 | i±1/2 and F 1 | i±1/2 , consistently on the basis of the CE solution of the Rykov model: for a well-resolved §ow, it is possible to obtain: where the derivatives are obtained analytically with the derivative of the Maxwellian de¦ned as follows: The time derivatives of the macroscopic variables can be obtained in terms of the space derivatives by means of the compatibility condition for the Rykov model: which, as explained in [40,41], leads to Introducing the CE solutions (7) of the Rykov model (1) with a modi¦ed collision time τ * : Since the di¨erence in the relaxation rate between translational and rotational processes is inherited in the collision number, Z r , let us de¦ne a single modi¦ed collision time by taking the moment relative to the total heat §ux.Thus, multiplying the ¦rst of Eqs.(10) by c ′ x c ′ 2 and the second one by c ′ x and adding the two resulting equations, one obtains: Let us simplify the numerical evaluation of this closure of the CE expansion by neglecting the terms relative to the Prandtl number correction introduced in the Rykov model; then, where ∂ ρT The second derivatives ∂ 2 /(∂t∂x) can be expressed in terms of only spatial derivatives thanks to the compatibility conditions ( 8) where W R and Q R are de¦ned in Eqs.(9), while the second time derivative of the mean velocity can be obtained considering the second-order compatibility conditions for the translational part: where as discussed in [40,41].
Finally, similar to [20,21], since DF M (T ) and D 2 F M (T ) will be sensitive to numerical errors (especially, close to equilibrium regions where they tend to vanish, and to impose that the physical stress inside the shock layer should be larger than the stress in the NS equation), a limiter is needed.In the current work, the following nonlinear limiter is used: where Kn is the local Knudsen number based on the gradients of the macroscopic variables de¦ned as in [12,42] and the function f (Kn) = −0.5 (1.0 + exp (−c Kn)) is used to obtain a smoother transition from τ * = τ to τ * = 2τ .Here, the parameter c = 4.01341 is calculated in order to have at least a starting τ * = 1.1τ at Kn = 0.05.

Di¨use Wall Boundary Conditions
The gas evolution at a solid boundary is modeled assuming that particles hit the wall with a distribution function according to the §ow conditions whereas they are re §ected with: a Maxwellian distribution according to the wall state for fully accommodation boundary (viscous wall); the same distribution function for specular re §ection boundary (inviscid wall); and a combination of di¨use and specular boundaries depending on the accommodation coe©cient σ.
Therefore, the ¦nal gas distribution function at the wall can be written as where u < 0 and u > 0 represent the velocities of particles hitting the wall and re §ected by the wall, respectively, while F M and F are the Maxwellian and the nonequilibrium distribution functions at the wall.The §uid state at the wall can be extrapolated from the domain.In the present work, fully accommodated walls have been employed, i. e., σ = 1.

RESULTS AND DISCUSSION
4.1 Assessment of the Gas-Kinetic Scheme Computational and memory cost is a major drawback for DVM in hybrid simulations and the use of the GKS to reduce the extension of the domain where  4)Eq.( 6) GKS it is strictly required represents a preferable alternative.To support this assertion, di¨erent test cases for a wide range of Mach numbers have been considered.Details of the simulations are reported in Table 1 and wall boundary conditions are considered fully accommodated.As shown in Figs.4a and 5a, in contrast to traditional single-/multitemperature NS approaches [44,45], the GKS is able to resolve shock structures with and without rotational nonequilibrium.However, due to the continuum formulation, it predicts steeper shocks relative to    DSMC solutions.Shock-structure predictions in better agreement with DSMC results can be obtained when Eq. ( 11) is employed to modify the collision time as shown in Fig. 5a.On the other hand, employing the modi¦ed collision time slightly a¨ects the accuracy in the subsonic region of the shock.The main differences between GKS and DSMC results can be observed for the higher Mach number cases in Figs.4b and 5b.This is probably due to the use of the CE expansion which does not allow to represent the typical bi-modal behavior of the distribution function across shock waves at very high Mach numbers (> 4).
The capability of the GKS to resolve multiple temperatures across shock structures, as shown also in Figs. 6 and 7 for the wedge and waverider cases, respectively, makes possible the reduction of the domain where more complex approaches, such as the DVM, are necessary.Indeed, as it is possible to observe from Figs. 811, the GKS results prove to be in good agreement with DVM and DSMC calculations with the biggest di¨erences between GKS and DVM, de¦ned here as occurring when the local Knudsen number (see Eq. ( 19)) is much higher than the commonly employed threshold, i. e., 0.05 [47], for the continuum breakdown.Furthermore, in both two-dimensional cases, the quantities at the wall are pre-   Velocity space cells size is equal to 0.5u∞.The DSMC, hybrid approach, and continuum results reported are from [11] dicted with a di¨erence less than 5% relative to DVM calculations and Fig. 11b also shows a good agreement between GKS and hybryd method results.Regarding the waverider case, as expected, the GKS predicts a §ow ¦eld with high thermal nonequilibrium around the waverider nose where also important viscous interaction e¨ects occur and in the boundary layer over the upper body surface (Figs. 7 and 12).For the cylinder test case, Fig. 13 shows density and temperatures pro-¦les for the section at 45 degree respective to the symmetry plane.The GKS correctly predicts the shock position as well as quantities and gradients at the wall in comparison to the results of the Modular Particle-Continuum (MPC) method of [42]; however, a thinner shock is predicted.A better shock structure prediction can be obtained using the hybrid approach as can be observed in Fig. 14.Here, for the hybrid approach, the DVM is employed in a region around the bow shock while the rest of the domain is simulated using the GKS (Fig. 15).The information exchange between the two solvers has  Mach 12 nitrogen §ow around a cylinder; velocity space cells size is equal to 0.5u∞.The MPC results (1) are from [42]; curves 2 refer to GKS based on Rykov model / DVM been handled employing the state-based coupling (see Fig. 3a), described in section 2 with an overlap region extension of about 10 cells.This means that across the overlap region, the distribution function employed as boundary condition for the DVM is reconstructed using the CE expansion (see Eq. ( 7)), where the macroscopic variables and the relative derivatives are de¦ned by the GKS solutions; vice versa: the macroscopic variables and the relative derivatives needed for the GKS §uxes are obtained from the DVM solution.The di¨erences can be still observed between the hybrid results from the present work and the ones in [42].The reason of these discrepancies can be found in the di¨erent methods employed.Indeed, in [42], a DSMC method is used where rarefaction e¨ect are dominant while in the present work, a DVM, which solves only a simpli¦ed model of the BTE where the collision time τ does not depend on the particle nitrogen §ow around a cylinder; velocity space cells size is equal to 0.5u∞ velocities, is employed.When employed in a hybrid simulation, as it is possible to notice from Fig. 16, the transition between the GKS and the DVM solvers at the interface is naturally smooth, no blending of the two solutions has been applied in the overlap due to the common root of the approaches.This and the extended validity of the the GKS also suggest a reduction of the hybrid simulations£ sensitivity to the positioning of the interfaces.However, the latter point requires further investigations.

Computational and Memory Cost
The GKS is implemented so that the particle velocity space dimensions depend on the local state, while due its complexity, this feature is not available in the DVM for the kinetic Boltzmann equations in MC.To perform a fair comparison of the computational time for the two approaches, constant velocity spaces were used.The runs have been performed on Intel R Xeon R processors at the University of Liverpool cluster ¤Chadwick¥ and solutions have been considered converged when the L 2 -norm of the update between two consecutive solutions is lower than 10 −7 for the normal shock cases and 10 −8 for the two-dimensional cases.Furthermore, it needs to be reminded that the halo exchange in the GKS involves only the §ow state while the DVM needs to exchange the full velocity space and this represents an advantage of the GKS relative to the DVM and the full UGKS when a parallel calculation is performed.In Table 2, the computa- tional time for the studied cases is reported and the GKS is found to be from 50% 90% faster than the DVM.This is partially due to the lower number of iterations needed in general by the GKS, but mainly to the time needed per iteration, around 710 times smaller than the DVM.Finally, in terms of the memory cost of the GKS, the latter is drastically reduced compared with the DVM and the full UGKS.Indeed, in the DVM and the UGKS, the values of the distribution function need to be stored for each physical cell in the full velocity space while the GKS being employed in the context of a continuum solver requires only the storage of the primitive variables.
Thus, employing the GKS in place of the DVM where the §uid is near thermal equilibrium, the performance of the hybrid solver can be improved in both memory and CPU time requirements, this without compromising the accuracy as shown in subsection 4.1.

CONCLUDING REMARKS
In the present work, an analytical GKS has been developed and added to the framework presented in [23,26,27].The prediction of §ow ¦elds where rare¦ed and continuum regions coexist requires the solution of two models; the NS equations and the BTE.Since the methods to solve the BTE are expensive, the reduction of the region where this is strictly required could improve the performance of hybrid simulations.For these reasons, a GKS for near-continuum regime has been proposed.
The scheme has been tested for various cases and Mach numbers proving to produce reliable predictions for near-continuum §ows.Moreover, the GKS also proved to be capable of solving more complex three-dimensional §ow ¦elds and to couple naturally with a DVM based on the same kinetic model.Regarding the computational time, when compared with a kinetic DVM solver, the nearcontinuum GKS solver was found to be between 50% and 90% faster than the former.Furthermore, due to the lower number of variables that need to be stored, the GKS is less expensive in terms of memory than the DVM and the full UGKS.This proves that GKS can be a viable way to improve the performance of hybrid simulations maintaining an acceptable level of reliability when in place of more complex methods for weakly rare¦ed §ows.

Figure 6
Figure 6 Mach 4 nitrogen §ow around a 25 degree wedge; velocity space cells size is equal to 0.5u∞: temperatures di¨erence, GKS solution

Figure 10
Figure 9The GKS and DVM results di¨erence (a) and local Kn (see Eq. (19)) (b).Mach 4.89 nitrogen §ow around a §at plate; velocity space cells size is equal to 0.5u∞

Figure 13
Figure 13  The ρ (a) and Tt (solid curves) and Tr (dashed curves) pro¦les (b) at 45 • .Mach 12 nitrogen §ow around a cylinder; velocity space cells size is equal to 0.5u∞.The MPC results (1) are from[42]; curves 2 refer to GKS based on Rykov model

Figure 15 Figure 16
Figure 15The DVM and GKS domains.Mach 12 nitrogen §ow around a cylinder; velocity space cells size is equal to 0.5u∞

Table 1
Test cases details

Table 2
Test cases details and computational time