# Towards Multidimensional Artificially Characteristic-Based Scheme for Incompressible Thermo-Fluid Problems.

Nomenclature

CFL - Courant-Friedrichs-Lewy, T - temperature, K, [C.sub.p]- pressure coefficient, u - velocity component in x direction, m/s, f(x,y,t) - characteristic surface, U- velocity vector (m/s), Fr - Froude number (Non-dim.), v - velocity component in y direction (m/s), [F.sub.c] - convective flux vector (Non-dim.), V - pseudo-velocity vector (m/s); [F.sub.v] - Viscous flux vector (Non-dim.), W - vector of dependent variables (Non-dim., Gr - Grashof number (Non-dim.).

Greek Symbols:

m - time level (Non-dim.), [beta] - artificial compressibility (Non-dim.), n - normal direction (Non-dim.), [[beta].sub.ex] - thermal volume expansion coefficient (1/K), Nu - Nusselt number (Non-dim.), [phi] - Generic name for flow variables (Nondim.), n - normal vector to characteristic surface (Non-dim.), [phi] - wave angle (rad.), [n.sub.t] - tangential-component of normal vector (Non-dim.), [eta] - connective parameter (Nondim.), p - pressure (Pa), [theta] - angular coordinate around cylinder (rad.), Pr - Prandtl number (Non-dim.), [rho] - density (kg/[m.sup.3]), q - heat flux (W/[m.sup.2]), [[tau].sub.ij] - Viscous stresses (N/[m.sup.2]), R - source term in Navier-Stokes equations (Non-dim.), [omega] - vorticity (1/s), Re--Reynolds number (Non-dim.), [psi] - defined as [psi] = [f.sub.t] +uf+v[f.sub.y] (Non-dim.), S - surface area ([m.sup.2]).

Superscript:

S' - aell area ([m.sup.2]), * - dimensional quantities (Non-dim.), t - tangential direction (Non-dim.).

1. Introduction

Various numerical techniques have been presented for incompressible flows with pressure correction schemes. The main difficulty is finding a relationship between velocity and pressure fields to satisfy the continuity equation. The artificial compressibility of Chorin can be considered as an alternative of these methods [1]. In this method, the continuity and momentum equations are coupled, thereby it is possible to use time-marching methods. Addition of artificial compressibility may affect the convergence process without interference on the results [2]. A few numerical schemes based on artificial compressibility have been developed by researchers. Among them, central scheme with Jameson's artificial viscosity was proposed to prevent the pressure-velocity decoupling [3]. The main idea of artificial compressibility methods such as Godunov-type schemes is similar to compressible methods. The Roe's flux difference splitting is used for incompressible flows [4-7]. They used artificial compressibility along with Roe averaging for flux discretization. Using pseudo-time derivative makes it possible to apply characteristic methods for incompressible equations. Conventional characteristic-based (CB) schemes which have been presented for fluid flows are constructed under the assumption of locally one-dimensional. Drikakis et al. extended this method for three-dimensional flows and incorporated it with multigrid techniques [8]. Zhao et al. utilized aforementioned method to unstructured two- and three-dimensional grids with heat transfer [9-13]. Flow variables are calculated along characteristic paths in the directions normal to the cell faces and their initial values are interpolated based on the sign of the corresponding characteristic speeds. The characteristic-based scheme has been used for a wide range of incompressible flows such as porous media [14], flow with various densities [15-16] and in conjunction with artificial compressibility [17-21]. The real multidimensional nature of flow demands to consider the directions in which the information is propagated. Different upwind schemes have been developed for the multidimensional Navier-Stokes equations. Razavi et al. introduced a genuinely multidimensional characteristic based scheme [22]. This scheme was developed for velocity field and tested for several benchmarks [23]. In the present study for the first time the two-dimensional characteristics for energy and Navier-Stokes equations are presented. The new multidimensional artificially characteristic-based (MACB) scheme is applied for square cavity with heat transfer and forced convection around a circular cylinder. The results were presented for a wide range of Reynolds and Grashof numbers. The main application of the proposed scheme would be in solving all the incompressible flows in different and complex geometries with heat transfer accurately such as heat exchangers, vortex tubes and etc.

2. Governing equations

The Navier-Stokes and energy equations for two-dimensional incompressible flow in non-dimensional form are expressed as:

[mathematical expression not reproducible] (1)

where

[mathematical expression not reproducible]

where [phi] is couples the energy, momentum and continuity equations. The viscous stresses and heat flux are given as follows:

[mathematical expression not reproducible] (2)

Based on the non-dimensional scaling, the above equations have been non-dimensionalized [9].

3. Mathematical formulation

To derive the characteristic relations for two-dimensional Navier-Stokes and energy equations, their equivalent Euler counterparts are considered. Characteristic equations can be obtained by assuming a characteristic surface as f(x, y, t) = 0.

[mathematical expression not reproducible] (3)

The kinematic relations relate the partial derivatives to total ones corresponding to the assumed surface and yield the following system of equations [24]:

[mathematical expression not reproducible] (4)

Where the subscripts stand for the partial differentiation.

To satisfy compatibility requirements of Equations (4), the determinant of coefficient matrix is set to zero, so

one has:

[mathematical expression not reproducible] (5)

Pseudo-velocity vector V = (u, v, 1) and normal to characteristic surface n = (cos [phi], sin [phi], [n.sub.t]), alike compressible Euler equations [25] are defined in which [phi] shows the wave direction. Equations (5) can be expressed in terms of V and n, thereby two types of characteristic surfaces corresponding to the following relations are obtained:

[mathematical expression not reproducible] (6)

The first equation of (6) expresses pseudo-streamlines and the second one demonstrates the virtual characteristic fronts. By some algebraic manipulations, [n.sub.t] takes the following form:

[mathematical expression not reproducible] (7)

The above roots appear to have always different signs. Also, the growth of influence and dependence zones around the pseudo-streamlines are revealed. Two roots of Eq. (7) demonstrate the biplanes tangent to characteristic surfaces passing through a certain point and producing Mach cones. By considering Eq. (7), it is concluded that dual characteristic surfaces would exist so that their corresponding tangential planes construct two Mach cones which extend from a certain point, namely domain of dependence and domain of influence.

Wave paths can be specified by an assumed characteristic path on f (x, y,t) = 0. Taking df/dt = 0 yields:

[mathematical expression not reproducible] (8)

Virtual acoustic waves correspond to the characteristic paths on the Mach cone. To obtain these paths, the above equations and the second relation of Eqs. (6) are combined to yield:

[mathematical expression not reproducible] (9)

Subtraction of above relations from each other yields:

[mathematical expression not reproducible] (10)

Hence, the generators of Mach cone from Eq. (10) can be expressed as:

[mathematical expression not reproducible] (11)

The compatibility equations are obtained by taking Eq. (4) and second value of [psi] in Eq. (5), as follows:

[mathematical expression not reproducible] (12)

It is possible to select different characteristic paths and their compatibility relations to estimate the cell interface values. The above-mentioned compatibility relations are used to model the convective fluxes, whereas viscous fluxes are calculated through averaging [26].

4. Convective fluxes

Fig. 1 demonstrates domain of influence which is the intersection of local Mach cone (or ellipse) with time level n. As seen in Fig. 2, four pseudo-waves are selected and Eqs. (13) are discretized along them to observe the physical behavior of domain of influence for face "*".

In Fig. 2, [[phi].sub.1] to [[phi].sub.4] are the angles between x -axis and normal vector to the presented directions. To evaluate convective flux at face "*", the ([[phi].sub.1], [[phi].sub.2]) pseudo-acoustic waves corresponding to n direction and ([[phi].sub.1], [[phi].sub.4]) waves corresponding to t direction are selected. There, n is normal vector to common face between (i, j) and (i+1, j) and t is normal to n. Different wave numbers and directions can be selected with corresponding compatibility relations for evaluating interface values. For choosing ([[phi].sub.1]) of face "*", the intersection of the line between midpoint of face "*" and the center of cell (i+1, j) with Mach ellipse is considered. The value of [[phi].sub.1] is determined by tangent to Mach ellipse at this point. This procedure is repeated for other waves. Using Eqs. (12), six compatibility equations corresponding to [[phi].sub.1] and [[phi].sub.2] are obtained as follows:

[mathematical expression not reproducible] (13)

Here, [n.sub.t1] and [n.sub.t2] obtained from Eq. (7) with respect to wave angles [[phi].sub.1] or [[phi].sub.2] where positive sign of [n.sub.t] is used. Similar compatibility equations can be written for waves [[phi].sub.3] and [[phi].sub.4]. Discretizing Eqs. (13) along with their corresponding paths, yields the following system of equations:

[mathematical expression not reproducible] (14)

where: A, B, C, D, E, F can be determined from Eqs. (7) and (13) as follows:

[mathematical expression not reproducible] (15)

and

[mathematical expression not reproducible]

Similar procedure is done for [[phi].sub.3] and [[phi].sub.4]. In Eqs. (14), p (*), u (*), ? (*) and T (*) denote the flow variables at cell interfaces. At first, u (*) and p (*)are calculated from the first and forth relations and the p (*) and ? (*) are determined using second and fifth relations. Then, T (*) and p (*)are calculated from third and sixth relations of Eqs. (14). An average of p (*) is then considered as final value. Each equation includes flow information which is transported from time m to time m+1 and none of them can be ignored. The calculated values for u (*) and ? (*) from Eq. (14) are projected onto direction n and similarly, the obtained values from corresponding [[phi].sub.3] and [[phi].sub.4] compatibility relations are projected onto direction t in order to find the velocity components on the cell interface ''*''.

4.1. Viscous fluxes and time integration

Viscous fluxes are computed by variables derivatives at the cell interfaces as shown in Fig. 3. For example, the first-order derivative at side AB is determined as follows:

[mathematical expression not reproducible] (16)

where S' is the area of AMBN. Typically, one may use

[mathematical expression not reproducible]

For time discretization, a 4th-order Runge-Kutta method was used.

The limitation of proposed scheme could be the numerical range of connective parameter. Also, viscosity has not been taken into account in deriving the characteristic equation. The benefits of applied research methodology include almost second order accuracy within lesser computed time.

5. Results and discussion

5.1. Square cavity with heat transfer

To validate the present MACB code, the results were compared with [27-29]. At solid walls, no-slip conditions are applied. The top and bottom walls are maintained at non-dimensional temperatures and side walls are thermally insulated. The u-velocity profile is presented and compared with CB and central scheme [28] benchmark solution in Fig. 4. The MACB scheme requires no artificial viscosity even at high Reynolds numbers. The results of MACB scheme for vorticity are compared with their counterparts where the agreement is remarkable Fig. 5.

The Bossinesq approximation is employed and Prandtl number was set to 0.71 in all cases. An appropriate range for [eta] is obtained by numerical tests and lies between one and two. To evaluate the accuracy and convergence of MACB scheme with conventional methods, a couple of simulations were conducted at different Reynolds and Grashof Numbers. Also, conventional CB was tested and its results were reported here. A typical error norm is defined as:

Fig. 6 demonstrates typical convergence histories of MACB and conventional CB scheme using the averaging method for energy equation at Re = 3162 and Gr = [10.sup.6] on the same grid structure. As seen, MACB scheme presents the most rapid convergence than averaging and conventional CB schemes.

[mathematical expression not reproducible] (17)

According to results, the maximum permissible CFL number for Re = 1000 and Gr = 100 on 80x80 uniform grid can be increased up to 1.6 for MACB, while it is limited to 0.9 for conventional CB.

The temperature fields and streamlines are depicted for Re = 1000 and Gr = 100 on 128x128 uniform grid and Re = 3162 and Gr = [10.sup.6] on 256x256 uniform grid in Figs. 7 and 8, respectively. The isotherms are clustered close to the bottom wall which points to the existence of steep temperature gradients in the vertical direction in this region.

These temperature gradients are weak on the other portion of cavity.

Fig. 9 shows the verification of the results for temperature profiles. In the limit, when the top wall is stationary, i.e. Re[right arrow]0, the corresponding temperature distribution approaches the linear profile obtained by conduction. In contrary, when the buoyancy effect is minor, Gr/Re << 1, in the middle portion of cavity the temperature changes are very small and rapid changes in temperature distribution occur in adjacent of top and bottom walls.

Figs 10 and 11 respectively refer to the results of u- and v-velocity profiles along horizontal line through the cavity centerline using MACB scheme in comparison with results of [31]. These figures demonstrate the change in basic flow character as Gr increases at a fixed Re.

5.2. Forced convection around a circular cylinder

In this section, steady and transient forced convection heat transfer around a circular cylinder was investigated using MACB as a benchmark. The far field boundary has intelligently been separated by the condition V. n which shows the flow direction. For flow around cylinder at

Re[less than or equal to]40, the regime is steady and increasing Re leads to appearing of transient nature of flow. Fig. 12 shows convergence history of MACB and convectional CB scheme using the averaging method for energy equation at Re = 100 on 80x80 grid. Fig. 1 demonstrates transient temperature contour at Re = 150.

Table compares the mean Nusselt number of MACB scheme and experimental results. Fig. 14 shows results of Comparison of the local Nu distribution along the cylinder surface between MACB and Baranyi [39]. Fig. 15 demonstrates pressure coefficient distribution along the cylinder surface.

6. Conclusions

Present study proposes a new multidimensional artificially characteristic-based (MACB) scheme for simulation of incompressible viscous flows with heat transfer. Multidimensional characteristic structure for energy propagation in incompressible flow is derived for the first time. The proposed MACB scheme is used in finite-volume form to evaluate convective fluxes with heat transfer in a square cavity for a wide range of Reynolds and Grashof numbers and forced convection around a circular cylinder. Also, for comparison purposes, the CB scheme with averaging for energy equation is used. It was found that MACB has remarkable faster convergence in comparison with CB scheme and averaging methods. At higher Richardson numbers, the conventional flux averaging was failed to converge properly. Obtained results using new proposed scheme are in good agreement with the benchmark solutions in the literature.

References

[1.] Chorin, A. J. 1972. A numerical method for solving incompressible viscous flow problems, Journal of Computational Physics 2: 12-26. http://dx.doi.org/10.1016/0021-9991(67)90037-X.

[2.] Kwak, D.; Kiris, C.; Kim CS. 2005. Computational challenges of viscous incompressible flows, Computers & Fluids 34: 283-299. http://dx.doi.org/10.1016/j.compfluid.2004.05.008.

[3.] Farmer, J.; Martinelli, L.; Jameson, A. 1994. Fast multigrid method for solving incompressible hydrodynamic problems with free surface, AIAA journal 32: 1175-1182. http://dx.doi.org/10.2514/3.12117.

[4.] Rogers, SE.; Kwak, D. 1991. Steady and unsteady solutions of the incompressible Navier-Stokes equations, AIAA Journal 29: 603-610. http://dx.doi.org/10.2514/3.10627.

[5.] Liu, C.; Zhang, X.; Sung, CH. 1998. Preconditioned multigrid methods for unsteady incompressible flows, Journal of Computational physics 139: 35-57. http://dx.doi.org/10.1006/jcph.1997.5859.

[6.] Kallinders, Y.; Ahn, HT. 2005. Incompressible Navier-Stokes method with general hybrid meshes, Journal of Computational physics 210: 75-108. http://dx.doi.org/10.1016/j.jcp.2005.04.002.

[7.] Yuan, L. 2002. Comparison of implicit multigrid schemes for three-dimensional incompressible flows, Journal of Computational Physics 177: 134-155. http://dx.doi.org/10.1006/jcph.2002.7007.

[8.] Drikakis, D.; Iliev, OP.; Vassileva, DP. 1998. A nonlinear multigrid method for the three-dimensional incompressible Navier-Stokes equations, Journal of Computational Physics 146: 301-321. http://dx.doi.org/10.1006/jcph.1998.6067.

[9.] Zhao, Y.; Zhang, B. 2000. A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids, Comput. Methods Appl. Mech. Engineering 190: 733-756. http://dx.doi.org/10.1016/S0045-7825(99)00443-0.

[10.] Tai, CH.; Zhao, Y. 2003. Parallel unsteady incompressible viscous flow computations using an unstructured multigrid method, Journal of Computational Physics 192: 277-311. http://dx.doi.org/10.1016/j.jcp.2003.07.005.

[11.] Tai, CH.; Zhao, Y.; Liew, KM. 2004. Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method, Computers and Structures 82: 2425-2436. http://dx.doi.org/10.1016/j.compstruc.2004.04.014.

[12.] Tai, CH.; Zhao, Y.; Liew, KM. 2005. Parallel-multigrid computation of unsteady incompressible viscous flows using a matrix-free implicit method and high-resolution characteristics-based scheme, Computer Methods in Applied Mechanics and Engineering 194: 3949-3983. http://dx.doi.org/10.1016/j.cma.2004.09.010.

[13.] Tai, CH.; Zhao, Y.; Liew, KM. 2005. Parallel computation of unsteady incompressible viscous flows around moving rigid using an immersed object method with overlapping grids, Journal of Computational Physics 207: 151-172. http://dx.doi.org/10.1016/j.jcp.2005.01.011.

[14.] Siong, K.; Zhao, CY. 2004. Numerical study of steady/unsteady flow and heat transfer in porous media using a characteristics-based matrix-free implicit FV method on unstructured grids, International Journal of Heat and Fluid Flow 25: 1015-1033. http://dx.doi.org/10.1016/j.ijheatfluidflow.2004.01.005.

[15.] Shapiro, E.; Drikakis, D. 2005. Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part I. Derivation of different formulations and constant density limit, Journal of Computational Physics 210: 584-607. http://dx.doi.org/10.1016/j.jcp.2005.05.001.

[16.] Shapiro, E.; Drikakis, D. 2005. Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part II. Multigrid implementation and numerical tests, Journal of Computational Physics 210: 608-631. http://dx.doi.org/10.1016/j.jcp.2005.05.002.

[17.] Drikakis, D.; Goldberg, U. 1998. Wall-distance-free turbulence models applied to incompressible flows, International Journal of Computational Fluid Dynamics 10: 241-253. http://dx.doi.org/10.1080/10618569808961688.

[18.] Mallinger, F.; Drikakis, D. 2002. Instability in three-dimensional unsteady stenotic flows, International Journal of Computational Heat Fluid Flow 23: 657-663. http://dx.doi.org/10.1016/S0142-727X(02)00161-3.

[19.] Drikakis, D.; Smolarkiewicz, PK. 2001. On spurious vertical structures, Journal of Computational physics 172: 309-325. http://dx.doi.org/10.1006/jcph.2001.6825.

[20.] Drikakis, D. 1997. Bifurcation phenomena in incompressible sudden expansion flows, Physics of Fluids 9: 76-87. http://dx.doi.org/10.1063/1.869174.

[21.] Neofytou, P.; Drikakis, D. 2003. Non-Newtonian flow instability in a channel with a sudden expansion, Journal of Non-Newtonian Fluid Mechanics 111: 127-150. http://dx.doi.org/10.1016/S0377-0257(03)00041-7.

[22.] Razavi, SE.; Zamzamian, K.; Farzadi, A. 2008. Genuinely multidimensional characteristic-basead scheme for incompressible flows, Int. J. Numer. Meth. Fluids 57: 929-949. http://dx.doi.org/10.1002/fld.1662.

[23.] Zamzamian, K.; Razavi, SE. 2008. Multidimensional upwinding for incompressible flows based on characteristics, J. of computational physics 227: 8699-8713. http://dx.doi.org/10.1016/j.jcp.2008.06.018.

[24.] Razavi, SE. 1995. Far field boundary conditions for computation of compressible aerodynamic flows, Ph.D. Thesis, Department of Mechanical Engineering, McGill University, Montreal, Canada.

[25.] Zacrow, MJ.; Hoffman, JD. 1976. Gas Dynamics, Vol. II, John Wiley and Sons, New York.

[26.] Ekaterinaris, JA. 2004. High-order accurate numerical solutions of incompressible flows with the artificial compressibility method, International Journal for Numerical Methods in Fluids 45: 1187-1207. http://dx.doi.org/10.1002/fld.727.

[27.] Erturk, E.; Gokcol, C. 2006. Fourth-order compact formulation of Navier-Stokes equations and driven cavity flow at high Reynolds numbers, Int. J. Numer. Methods Fluids 50: 421-36. http://dx.doi.org/10.1002/fld.1061.

[28.] Ghia, U.; Ghia, KN.; Shin, CT. 1982. High-resolutions for incompressible Navier-Stokes equation and a multi-grid method, J. of Computational Physics 48: 387-411. http://dx.doi.org/10.1016/0021-9991(82)90058-4.

[29.] Bruneau, CH.; Saad, M. 2005. The 2D lid-driven cavity problem revisited, Comput. Fluids 35: 326-48. http://dx.doi.org/10.1016/j.compfluid.2004.12.004.

[30.] Botella, O.; Peyret, R. 1998. Benchmark spectral results on the lid-driven cavity flow, Computers and Fluids 27: 421-433. http://dx.doi.org/10.1016/S0045-7930(98)00002-4.

[31.] Iwatsu, R.; Hyun, JM.; Kuwahara, K. 1993. Mixed convection in a driven cavity with a stable vertical temperature gradient, Int. J. Heat Mass Transfer 36: 1601-1608. http://dx.doi.org/10.1016/S0017-9310(05)80069-9.

[32.] Cheng, TS.; Liu, WH. 2010. Effect of temperature gradient orientation on the characteristics of mixed convection flow in a lid-driven square cavity, Computers and Fluids 39: 965-978. http://dx.doi.org/10.1016/j.compfluid.2010.01.009.

[33.] Karamers HA. 1946. Heat transfer from spheres to following media, Physica 12: 61-80. http://dx.doi.org/10.1016/S0031-8914(46)80024-7.

[34.] Fand RM. 1965. Heat transfer by forced convection from a cylinder to water in crossflow, Int. J. Heat Mass Transfer 8: 955-1010. http://dx.doi.org/10.1016/0017-9310(65)90084-0.

[35.] Whitaker S. 1972. Forced convection heat transfer calculations for flow in pipes, past flat plate, single cylinder, and for flow in packed beds and tube bundles, AIChE J. 18: 361-371. http://dx.doi.org/10.1002/aic.690180219.

[36.] Perkins HC.; Leppert G. 1962. Forced convection heat transfer from a uniformly heated cylinder, J. Heat Transfer 84: 257-263. http://dx.doi.org/10.1115/1.3684359.

[37.] Perkins HC.; Leppert G. 1964. Local heat-transfer coefficients on a uniformly heated cylinder, Int. J. Heat Mass Transfer 7: 143-158. http://dx.doi.org/10.1016/0017-9310(64)90079-1.

[38.] Churchill SW.; Bernstein M. 1977. A correlating equation for forced convection from gases and liquids to a circular cylinder in cross flow, J. Heat Transfer 99: 300-306. http://dx.doi.org/10.1115/1.3450685.

[39.] Baranyi L. 2003. Computation of unsteady momentum and heat transfer from a fixed circular cylinder in laminar flow, Journal of Computational and Applied Mechanics 4: 13-25.

A. ROSTAMZADEH (*), S. E. RAZAVI (**), S. M. MIRSAJEDI (***)

(*) Department of Mechanical and Aerospace Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran

(**) School of Mechanical Engineering, University of Tabriz, Tabriz, Iran, E-mail: razavi@tabrizu.ac.ir

(***) Aerospace Eng. Dept., Faculty of New Technologies & Engineering, Shahid Beheshti University, G C., Tehran, Iran

crossref http://dx.doi.org/10.5755/j01.mech.23.6.15804

Received July 28, 2016

Accepted December 07, 2017
```Table
Comparison of mean Nusselt number for laminar flow passing through a
circular cylinder

Range of  Author                    Mean Nu
Re                                  number

Karamers [33]             2.00
Re = 10   Fand [34]                 1.79
Whitaker [35]             1.34
MACB                      1.72
Karamers [33]             2.67
Re = 20   Fand [34]                 2.45
Whitaker [35]             1.93
MACB                      2.35
Karamers [33]             3.62
Perkins and Leppert [36]  2.67
Re = 40   Fand [34]                 3.39
Perkins and Leppert [37]  2.83
Whitaker [35]             2.80
MACB                      3.25
Karamers [33]             3.99
Perkins and Leppert [36]  3.03
Re = 50   Fand [34]                 3.78
Perkins and Leppert [37]  3.21
Whitaker [35]             3.16
MACB                      3.63
Karamers [33]             6.64
Perkins and Leppert [36]  5.67
Fand [34]                 6.52
Re = 150  Perkins and Leppert [37]  6.03
Whitaker [35]             5.72
Churchill and Bernstein   6.22
MACB                      6.21
```
COPYRIGHT 2017 Kauno Technologijos Universitetas
No portion of this article can be reproduced without the express written permission from the copyright holder.