Fuzzy prediction of AWJ turbulence characteristics by using typical multi-phase flow models

ABSTRACT For the purpose of improving the working performance of AWJ (Abrasive waterjet) grinding in actual machining, and revealing the unknown influence mechanism of waterjet machining prediction, this research conducted a theoretical and experimental investigation concerning with the influences caused by multi- phase flow models, on the fuzzy prediction of turbulence characteristics. First, typical flow models describing the abrasive waterjet were presented and enunciated; Second, a series of turbulence characteristics were defined to quantitatively demonstrate the waterjet flow, which are responsible for the actual performances of AWJ grinding in the contact zone; Then a new fuzzy prediction system equipped with logic protocols was established and introduced to predict turbulence characteristics based upon the machining parameters by employing different flow models, therefore a set of quantified prediction results generated from the CFD (Computational Fluid Dynamics) simulation, orthogonal experiments and actual measurements can be obtained and distinguished systematically. Through a detailed performance analysis with benchmark indexes and prediction comparisons, the theoretical superiorities and specific applications of these flow models were fully discussed and then an integrated conclusion assessing their inherent influence mechanism was acquired. Based on these innovative achievements, the instantaneous machining inspection or parameter optimization in AWJ grinding process can be facilitated, simultaneously the new research ideas aimed at waterjet monitoring or grinding improvements can be presented as well.


Introduction
For industrial applications, AWJ grinding was widely used for workpiece profiling or surface machining. It is noteworthy that abrasive waterjet flow exerts remarkable influences on the grinding effectiveness of objective workpieces and the mechanical properties of material specimens in the contact zone. During the process of AWJ grinding, abrasive grits were combined with waterjet slurry therefore the particle-laden flow make a great contribution to the topography transformation on targeted surface. But as the influence mechanism caused by multi-phase waterjet flow on the turbulence characteristic determination still keeps in an obscure situation, this article decided to make a theoretical and experimental study concentrating on this topic.
Through searching the literatures published in the domain of AWJ machining, Bilbao Guillerna, Axinte, and Billingham (2015) have proposed some innovative ideas on the energy stream of abrasive waterjet; meanwhile, a directly numerical simulation, current AWJ researches were regrettably weakened in accuracy and reliability.
On the topic of the performance analysis between different flow models, Becker (2013a, 2013b); Li et al. (2013) respectively focused on the finite element modeling of waterjet flow, and the radial-mode AWJ machining process as well. More progresses relating to the important influences caused by flow models on the AWJ characteristic determinations can be found from (Atashbar Orang, Razavi, & Pourmirzaagha, 2014;Cui, An, Mao, & Li, 2009;Kumar Pal & Choudhury, 2014;Lebar, Junkar, Poredoš, Cvjeticanin, 2010). After introducing the idea of fuzzy prediction, we can see that more researchers began paying their high attentions on this new method, since it returns an accurate and stable prediction result from miscellaneous data sets. For example, Erkan, Demetgül, Işik, and Tansel (2014) used the Taguchi method and GONNs (Genetically Optimized Neural Network systems) to analyze the discrete data of abrasive flow and the fuzzy selection of waterjet parameters; Zain, Haron, and Sharif (2011) employed the ANN (Artificial Neural Network) and annealing approaches to predict the optimal process parameters of flow models. Ergur and Oysal (2015) studied on the fuzzy optimization of waterjet machining, and the estimation of cutting speed in abrasive waterjet either. Furthermore, Shabgard, Badamchizadeh, Ranjbary, and Amini (2013) investigated the fuzzy predictions of material removal rate, tool wear ratio, and surface roughness, with typical flow models were carefully considered about. Labib, Keasberry, Atkinson, and Frost (2011) have proposed a fuzzy logic controller for electrochemical machining (ECM). Meanwhile, literatures (Franco, Gasche, & Neto, 2016;Liang, Xie, Liao, & Zhou, 2015;Moussa, Della Valle, Garnier, & Peerhossaini, 2016;Pan, Ko, Yang, & Hsu, 2011) also reported their latest progresses on some relevant methods and experimental applications. But these investigations hardly covered all involved participant elements in our research, ranged from the appropriate turbulence characteristics for waterjet description, to the multi-phase flow models employed in grinding evaluations. There was still no literature discussed the fuzzy influence prediction of multi-phase waterjet flow, therefore the fundamental problems raised earlier remain unaddressed. Fortunately, it is worth noting that Chau et. al. have made great contributions to the actual applications of contemporary soft computing techniques in theoretical fields, including the predictions of daily flows, river stage, and plant leaf classification, etc. (Bhaumik, Bhattacharyya, Nath, & Chakraborty, 2016;Wang, Chau, & Xu, 2015;Wu, & Chau, 2011;Wu, Chau, & Li, 2009;) Therefore a series of optimal prediction results characterized by higher precision and better stability were obtained, which provided an important theoretical basis and reliable tool for the establishment of fuzzy prediction system and soft computing techniques undoubtedly; Moreover, based on the enlightenment of contemporary soft computing, the mechanism investigation of turbulence characteristic predictions can be greatly improved, too (Taormina & Chau, 2015). As the fuzzy prediction of AWJ characteristic streamlines the machining efficiency and controllable monitoring in industrial practices, it should be discussed thoroughly.
In this article, Section 1 explained the objective of the fuzzy prediction of AWJ turbulence characteristics by using typical multi-phase flow models, and the investigations of influence mechanisms caused by the latter; Section 2 gave the governing equations of typical flow models; Section 3 presented the representative mathematical characteristics to describe abrasive waterjet turbulence; and then Section 4 proposed the details of dynamic CFD simulations and actual AWJ grinding, Section 5 demonstrated the fuzzy prediction system being used. Thereafter an orthogonal experiment verification and influence mechanism evaluation for the fuzzy prediction of turbulence characteristics were completed in Section 6. Finally, Section 7 concluded this article as required.

Governing equations of typical multi-phase flow models
Since turbulence flow is very difficult to be measured directly, typical multi-phase flow models should be used for flow demonstration and principle investigation. For the convenience of this research, we aimed at studying their influences on the prediction of turbulence characteristics, with k-, k--R t , Spalart-Allmaras, Reynolds stress and large eddy simulation (LES) models were focused sequentially.

k-model
The standard k-model is written in each phase p as: With C ε1 = 1.44, C ε2 = 1.92, C ε4 = 1.2, p = f/s denotes the fluid or solid phase. D μ εp is the kth diffusion term : The fluid viscosity is μ T p = ρ p C μ (k 2 p /ε p ) with C μ = 0.09. Based on this arrangement the waterjet turbulence in dispersed phases was directly related to that of continuous phase, through the turbulence response coefficient p μ p is the production due to the mean velocity gradients (Bankwitz, Bankwitz, Förste, & Schulze, 1990): Here v p represents the velocity of fluid or solid phase, α p denotes the volume fraction of each phase, p i p is the specific two-phase turbulence term, represents the influence of the ith phase on the pth phase, p b p is the single phase buoyancy term modeling the gravity and density gradients on fluid or solid phase. In the th equation, max (p b p , 0) stands for a stable stratification if p b p ≤ 0; or for an unstable stratification if p b p ≥ 0. σ t is the turbulent Prandtl number equal to 0.9 by default (Delele, Weigler, Franke, & Mellmann, 2016).

k--R t model
k--R t model consists by the following expressions and equations: Here, μ p denotes an undamped pseudo-eddy viscosity, k p , c kp , n p , δ ij are specific coefficients of fluid or solid phases according to the fluid Reynolds stress equations.
expresses the fluid Reynolds stress production/destruction: The physical meanings of k--R t models are (a): the particle-laden turbulence fluctuation is determined by convection, diffusion and production; (b): the particleladen turbulence fluctuation is anisotropic and stronger than that of fluid turbulence. Through experimental comparisons it can be seen that k--R t model gives a more uniform particle concentration and spatial distribution (Zhou, 2010).

Spalart-Allmaras model
In Spalart-Allmaras turbulence model, the transport equation for viscosity quantityμ is, uses this quantity to modulate fluid viscosity. ρ f represents the flow density, μ m = μ t /ρ f denotes the molecular kinematic viscosity, and μ t denotes the molecular dynamic viscosity (Adel, Akhtari, & Dehghani, 2016). Its transport equations are provided below: Here s is defined as s gives the magnitude of vorticity. The constant d is the distance from the targeted cell to the closest boundary wall. On the walls,μ = 0.
As the nominal kinematic viscosity is constant (Das et al., 2015), it is modified tõ The spatial derivatives ofμ t , v t x and v t y can be defined as: ∂ 2 s ∂y 2 = s y+ y − 2s y + s y − y y 2

Reynolds stress model
In this section Reynolds stress model includes the waterjet of two-directional coupling between the solid particles and fluid phases, and the transportation of incompressible fluid (the carrier fluid) laden with abrasive particles (the solid phase) was considered about. The volumeaveraged particle equation can be described as: For the flow-phase Reynolds stress: For the solid-phase Reynolds stress: The production term is calculated as The source term R f accounts for the fluid-induced turbulence, the total turbulence kinetic energy is calculated from the trace of the total Reynolds stress: (Jovanović, Miljanović, & Jovanović, 2015). The Reynolds stress is calculated as following: Here C vm = 1.2, a = 0.1 and b = 0.3. For the fluidparticle velocity covariance.
is the pressure-dispersed phase velocity gradient correlation. To close Reynolds stress model, the transport equation (Nimvari, Maerefat, Jouybari, & El-hossaini, 2013) of the dissipation rate f should be solved and expressed as: Here the source term S εf is set to zero in this work. The model constants are listed as C ε1 = 1.44, C ε2 = 1.92, C s = 0.15 (Liu, & Hinrichsen, 2014;Taulbee, Mashayek, & Barre, 1999).

Large eddy simulation (LES)
LES was frequently used to validate the two-phase flow model, both fluid (coolant fluid or jet water) and solid phases (grinding grits or abrasive particles) were treated as interpenetrating continua. The filtered governing equations for LES model are given: For particle-collision stress: The definition of fluid viscosity μ t is ), which describes the anisotropic flow problems. Here is the length scale (or grid filter) defined in terms of the cell volume V, locally changed as = min(k, C w V 1/3 ); k is the von Karman constant, d is cell distance from the wall and C w = 0.544 is the switch coefficient. The strain tensors are defined as The energy equation is expressed as: The fluid Prandtl number was set as 0.9, which corresponds to a heat conductivity coefficient of 0.0257 W/m.K (Zhou, 2010). The momentum equation for the solid phase p (p = 1,2, . . . ,M) can be rewritten as: Here the terms on right-sides stand for the gravitational force without buoyancy effects, the pressure gradient, the molecular viscous, the momentum exchange between fluid and solid phases due to drag force, and moment exchange between them. pm is the particle-fluid momentum transfer coefficient, which may be recognized as the momentum source. It is assumed to be a function of particle fluctuation velocity.
The particle stress tensor is expressed in Newtonian form as (26) In this research, semi-empirical correlation was adopted based on the experimental determination of an effective viscosity coefficient for particulate flows, where the viscosity is proportional to the volumetric concentration of solid particles; v s is solid particle velocity, I is unit tensor, μ s = c s s and μ b = 0. Furthermore, for the fluid phase, Newtonian stress tensor is appropriate to describe viscous dissipation for the fluid phase: Here, v f is fluid velocity, the Reynolds stress τ ft can be described by LES approach through the Smagorinsky closure (Mason, 1994): Here, μ ft = ε f ρ f l 2 |S| is called as the fluid viscosity and S is the resolved deviatoric stress. The filter length l is proportional to the grid size away from boundaries l = C s , Cs = 0.1.

Abrasive turbulence characteristics
Since the high-speed particle-laden waterjet is very difficult to be identified by ocular inspection, a quantitatively demonstration of turbulence characteristics was required, which offers an important tool to monitor the AWJ grinding accurately. As some frequently used flow characteristics, such as the turbulence RMS (Root Mean Square) velocity (m/s), turbulence intensity, turbulence kinetic energy (KJ), turbulence entropy (J/K), and Reynolds shear stress (KPa), remarkably relate to the machining quality and grinding efficiency on specimen surfaces, they were investigated in details by this research.

Turbulence RMS velocity (V RMS )
According to the Newton laws of motion, the equation of motion for solid particles can be written as Here v f is the fluid phase velocity, v s is the solid particle velocity, μ f is the molecular viscosity of fluid, ρ s is the solid particle density, d s is the particle diameter, R e = (ρ s d s |μ s − μ f |/μ f )is the relative Reynolds number, and C D is drag coefficient. The instantaneous velocity of surrounding fluid is u f . The particle drag coefficient is: The grit movement is dictated for position and velocity, its position can be quantified by (d x s /dt) = v s Here x s = (x s ,y s ,z s ) is the solid article location and v = (v x ,v y ,v z ) is the instantaneous particle velocity vector. The RMS velocity of turbulent flow can be determined with a spherical particle of diameter d s (Vendra, Wen, & Tam, 2013): Here g is the acceleration due to gravity. The response time: Since this controls the particle dynamics, the time step independence is ensured by employing t < < T.

Turbulence intensity (T i )
Here the solid particles obey the local follow mode, T s stands for the intensity of solid particle pulsation, T s = (1/3)v 2 s , and v s denotes the RMS velocity of one given abrasive particle. Then the governing equation of multiphase turbulence can be given as (Terashima, Sakai, Nagata, Ito, & Onishi, 2014): Here α k , ρ k denote the volume concentration and density of phase k, k = f/s stands for the fluid phase or solid phase respectively, J is the Jacob determinant, g fm is the contravariant metric tensor; U k f is the contravariant component of phase k; υ k is the diffusion factor of k , and which represents U k i (i = 1,2,3), K f , f and T s ; U k i is the covariant RMS velocity component of phase k, K f is the kinetic energy of fluid phase; f is the dissipation rate of fluid-phase kinetic energy; T s is the temperature of solid phase, and S k stands for the source term of K f equation, it can be shown that: Here G k is the generate term of turbulence kinetic energy, D e is the dissipation term of fluid-phase turbulence intensity caused by the connection between abrasive particles and fluid flow. For the equation of f , its source term S f can be expressed as: Here C 1 and C 2 are empirical constants. The increment of turbulence intensity resulted from the solid particle wakes can be identified as: Here is empirical constant, v f and v s are the RMS velocity of fluid and solid phases. In the case of the density ration between solid particle and fluid flow is much large, K f is mainly determined by gravity force (Di Martino, Loia, & Sessa, 2011), thus: Here d s stands for the diameter of abrasive particle, and g denotes the gravitational acceleration. Then: Since the existence of abrasive particle changes the velocity distribution of fluid phase on the boundary layer a lot, the coefficients of molecular diffusion are considered to solve the boundary function: Based on this precondition, the turbulence intensity of multi-phase waterjet can be described as: Herev f ·v s denote the mean RMS velocities of fluid and solid phases, and T i represents the turbulence intensity.
gives the standard deviation, M and N stand for the numbers of experimental samples for fluid and solid phases. P(v) is the mean value of waterjet energy power (Vendra et al., 2013).

Turbulence kinetic energy (T c )
The fluid-phase motion is formulated in a manner similar to that of the usual model, except for the inclusion of particle concentration and the interaction forces between fluid and particles. The transport equation of dissipation rate for two-phase turbulence kinetic energy is: Here the new source term is Gs , D p,ij , P p,ij , p,ij are the diffusion, production terms of Reynolds stress and the production term concerning with fluid turbulence, respectively. For a closed system, the transport equations of n s v si , n s v sj , n s n s , v si v j , v sj v i were used. The fluid and particle momentum equations were closed as: The last term on the right side of the above-mentioned equations is closed, as the dissipation of velocity correlation keeps proportional to that of turbulence kinetic energy.
All these equations constitute the unified second-order moment of two-phase turbulence model. The balance for the fluid turbulence kinetic energy: Here v fi is the fluid fluctuating velocity based on k-model, and σ k , μ ef , μ f and μ s are the effective fluid viscosity, eddy viscosity, flow Prandtl number, and particle-phase viscosity, respectively (Kuwata, Suga, & Sakurai, 2014).

Turbulence entropy (T e )
Entropy has been applied to clarify the dynamic characteristics of flow patterns, it can be used as an efficient diagnostic tool to observe the flow transition. Therefore, this characteristic was introduced to identify the flow pattern and reveal the dynamic characteristics of flow pattern from macro and local views.
Since the irreversible physical process in energy conversion must be accompanied by energy dissipation, turbulence entropy has its own transport equation for singlephase incompressible flow (Franco et al., 2016).
Here s is the specific turbulence entropy; v x , v y and v z are the velocity components along the x, y and z directions in Cartesian coordinate system, respectively;q is the heat flux density vector; and T is temperature. In the Reynolds time-averaged process, the instantaneous variable can be separated into two parts,s = s + s , namely, the mean quantity part and the fluctuating part. Thus, the transport equation can be changed as follows.
For the incompressible flow, the viscous dissipation function is: As for turbulent flow, however, the specific entropy production rate can be separated into two terms after the Reynolds time-averaged process: one with a mean term and the other with a fluctuating term S pro,D = S pro,VD + S pro,TD Here S pro,VD is the specific entropy production rate caused by the time-averaged movement, called as the directly entropy production; S pro,TD is the specific entropy production rate caused by the velocity fluctuation, called as the turbulence entropy production. Both of them represent the local entropy production rates and can be defined as follows: The local entropy production terms can be calculated directly with 3D flow simulations by using ANSYS, and the overall entropy production of the computational domain was calculated by volume integration. After a detailed analysis of entropy production, it can be calculated directly when the 3D flow field was obtained. The entropy production theory has its own advantages for evaluating the rheology performance of abrasive waterjet grinding. Furthermore, this characteristic offers a quantitative and reliable analysis of turbulence energy dissipation for waterjet monitoring.

Reynolds shear stress (R s )
In this research, fuzzy prediction was based on the steady state solutions of Reynolds-averaged Navier-Stokes equations, which expresses the conservation of mass and momentum. The Reynolds shear stress can be obtained directly from: Here x denotes the axial distance, y denotes the lateral distance, and z denotes the span-wise one. δ ij is Kronecker delta, and represents the dissipation rate of turbulence kinetic energy. In the above equation, k is von karman constant, A ij stands for the term of re-distributive fluctuating pressure, while A ij w is the correction term for wall reflection effects. In which i, j, k, l,m are the indexes of coordinate axes, V, v' present the mean and fluctuating velocities in x-direction.
With D/Dt = ∂/∂t + V j ∂/∂x j , the term of redistributive fluctuating pressure is given by: Here A ij is modeled as a linear function of Reynolds shear stress tensor. The model constants were taken as Cs = 0.18, C 1 = 3.00, C 2 = −0.44, C 3 = 0.46, and C 4 = −0.23 (Rhea, Bini, Fairweather, & Jones, 2009). The turbulence energy dissipation rate is obtained from: With the standard constant values of C = 0.15, C 1 = 1.44 and C 2 = 1.90. Figure 1 presents the waterjet grinding machine employed in our research. In the condition of AWJ grinding, the machining performance depends on the interrelationships between different influence factors, especially those multi-phase turbulence characteristics should be paid high attention to. Since AWJ grinding can be improved considerably by the monitoring of turbulence characteristics, and the establishment of knowledge system with flow models, it is important to know the universe boundaries between them in advance. In this research, the actual AWJ grinding and CFD simulations were combined together to deal with this problem. Figure 2 gives the structure sketches of AWJ grinding head and monitoring devices. The AWJ grinding operation was undertaken by a high-pressure intensifier pump and 5-axis freeform grinding machine, (Zeeko IRP200 from ZeekoTM Ltd. of UK), with its construction details or parameter setting can be learned from the published literatures of Liang, Liu, Ye, and Wang (2013); Liang et al. (2015); Pang, Nguyen, Fan, and Wang, (2012). The outer rings of GCr15 deep groove ball bearing (DGBB) by the diameter of 72.00 mm were used for tests. This experiment controlled the particle concentration in the whole grinding process, and the abrasive particles uniformly distributed within the slurry flow. As the selected nozzle has the inner diameter of 0.50 mm and chamber length of 75 mm, only those major and controllable variables were studied, including the Waterjet pressure (P w /MPa),   n this research, Computational fluid dynamics (CFD) methodology was applied to study the influence mechanism of flow models. In order to demonstrate the waterjet flow and its machining process, a 3D numerical model has been set up to show its dynamic variation. Volume Of Fluid (VOF) was frequently used since it is particularly suitable to track the flow paths instantaneously. Since VOF simulation strongly depends on the grid resolution, a grid refinement should be prepared in the studied domains (i.e. The abrasive waterjet flow and waterjet grinding space), to improve its accuracy and decrease the computational effort (Chen, Zhang, & Feng, 2016). It is noteworthy to point out that the VOF method is an interface-resolving method and has become a popular methodology for solving two-phase flow phenomena. It was implemented to determine the interface between the fluid and solid phases, then the interface between different phases can be tracked by the volume fraction of fluid phases at the computational cell volume schemes. Based on this presupposed condition, water downstream pressure was arranged from 200 to 300 MPa for machining, while atmospheric pressure was presupposed as the boundary condition. The no-slip-shear was set as the boundary condition of flow wall, which means that fluid in contact with the machined surface got zero velocity (Zhuo & Zhong, 2013). The attack angle of abrasive waterjet varies according to the given workpiece specimen.

CFD simulations and actual AWJ grinding
Since the resolution of computational mesh impacts on the validation of numerical results greatly, the gridconvergence was outlined to choose an optimal grid for the fuzzy prediction and estimate the possible discretization errors (Gholami, Bonakdari, Zaji, & Akhtari, 2015). The simulations were run to generate a steady waterjet profile with the same experimental conditions. A set of 3-D hybrid mesh grids were investigated initially, the computational domain at each grid scheme keeps an identical proportion to the grid number. According to the actual experiment, the mesh cell sizes were tested and investigated by the sequence of 0.15,0.20,0.25,0.30,0.35,0.40 and 0.45 mm. It is also important to point out that, all these grids were obtained by strict refinements. All computations were initialized by using free-stream values (Isenmann et al., 2016). The grid number, computation efficiency and computation error were selected out as the verification indexes, with Table 3 shows the comparison results. For the sake of stability and accuracy, the grid scheme 3, with the unit size of 0.25 mm, was eventually selected out for the numerical calculations and machining description, because it returns a more reliable  (Liang, Zhou, Liu, & Wang, 2014), and the rest of simulation settings was reported by Table 4. In view of the obtained results, a satisfactory grid keeps a reasonable balance between accuracy, calculation efficiency and computational cost. Based on this setting preparation, the obtained mesh grid will be employed in sequence, and then Figure 3 describes the establishment of model grids and the turbulence diffusion on workpiece surface. Figure 4 demonstrates the AWJ grinding simulation, with the targeted areas in the contact zone of workpiece surface were highlighted by black lines, which describes the turbulence characteristics and its dynamic variation before waterjet flow touches the specimen surface. Based on the setup arrangement and computer simulation, Figure 5 gives the diffusion process of waterjet turbulence on the objective ring surface, from the contacting, machining, to the diffusing stage. From the  detailed CFD dynamic simulation of AWJ grinding, the planar distributions of multi-phase turbulence characteristics on the targeted contact areas were obtained and illustrated by Figures 6-10. The positions of turbulence characteristics in unstable situation were highlighted, with the same as the concentration regions of abrasive particles too (Errico, & Stalio, 2014). As it was difficult to quantify them by ordinary methods, we introduced a polar system of coordinates f(r,θ ) as Figure 11 illustrates. The radial coordinate is denoted by r, and the angular coordinate by θ . f(r,θ ) presents the turbulence characteristic values at the specific place of targeted coordinate point (Nakkina, Arul Prakash, & Saravana Kumar, 2016). Thereafter the averaged value of turbulence characteristics T characteristic can be calculated as: As an exampled result, Tables 5-9 demonstrate the obtained T characteristic values in some representative cases, thereafter their value distributons can be quantified correspondingly.

Fuzzy prediction
To establish an effective fuzzy prediction system for determining the AWJ turbulence characteristics with typical flow models, the value universes of input and output variables should be portioned accordingly (Londhe, 2008). Figure 12 demonstrates the detailed arrangements of two exampled membership functions, as the waterjet pressure (input) was illustrated by Figure 12(a) and partitioned from 0 to 8, or the turbulence intensity (output) in Figure 12(b), in which any value beyond the presupposed range is infinity or zero. 0 was assumed as 'Very Low' and 8 as 'Very High,' other value levels can be identified by 'Low,' 'Medium,' and 'High' in sequence (Liang, Liu, Tan, & Wen, 2016).
Then the fuzzy logic protocols were defined as: and "x 2 " is "A i2 " and "x 3 " is "A i3 " and . . . "x r "is "A ir ", then "f i " is "B i " (for i = 1, 2, 3, . . . , k n ) Here f' denotes the predicted output, p 1 , p 2 , . . . ,p r and r i are the consequent characteristics of the ith logic protocol. A i1 , A i2 , A i3 , . . . ,A ir are the linguistic labels representing the fuzzy data sets, whose membership functions are premise characteristics. x 1 , x 2 , x 3 . . . and x r denote the input variables; and k n represents the number of logic protocols. Here i , j , k , l and m were the predicted value levels; this expression represents the protocol antecedent, Figure 13 illustrates the exampled logic operation of Protocol 1 Figure 5. The diffusing simulation of waterjet turbulence on the objective ring surface.
and Protocol i. The output membership for multi-phase turbulence characteristics were matched and truncated according to fuzzy logic protocols: Here l denotes the selected membership function (from very high to very low), all truncated membership functions (ψ TurbulenceCharacteristics ) were aggregated by superposing the triangular shapes to create the final function (ψ out ) (Amano, Arakawa, & Suga, 2014): Table 10 demonstrates some representative fuzzy logic protocols by using Spalart-Allmaras model, Figure 13 shows the detailed aggregation in graphical illustration, and Figure 14 gives the predetermined protocols enunciated graphically.
The next step is defuzzification, in which the turbulence characteristics were computed by the centroid      method, it returns the center of shaped area in illustration figures under the curves and lines according to Eq. (54), with which the predicted turbulence characteristics can be obtained (Mohamad, Zain, Bazin, & Udin, 2015): Based on the obtained index of T C (i ), it demonstrates the predicted characteristics in the portioned universe ranged from 0 to 8, thereafter the turbulence characteristics can be transferred into the predicted value (Liang, Liu, & Ye, 2012):   Here i = 1,2,3,4,5, which represent V RMS , T i , T c , T e , and R s respectively. In this work, the fuzzy prediction of turbulence characteristic depended upon the past observations by historical measurements, which provided better predictions and more reliable results over other traditional works.

Orthogonal experiment for the verification of fuzzy prediction
In order to evaluate the influences of multi-phase flow models on the turbulence characteristic prediction, a series of machining tests were operated. Figure  15 exampled the performance comparisons of abrasive  waterjet grinding on the ring surface, by which a remarkable surface quality improvement can be achieved. In order to prove the accuracy and reliability of these predicted results, an actual measurement of waterjet variation was conducted by a Dantec Laser Doppler Anemometry (LDA) system, also known as the Laser Doppler Velocimetry (LDV). In this experiment, the high-resolution LDA monitoring device was mounted at the right edge of machining space, and employed to track the dynamic motion and disequilibrium propagation of waterjet flow on the machined surface, with the laser scanning by 20-30 Frames Per Second (FPS) simultaneously. After that, a vector diagram of waterjet flow and a spatial representation for fluid bouncing distribution can be acquired, which contributed to the integrated expressions of dynamic waterjet variations and energetic changes. Since the real-time measurements of multi-phase abrasive waterjet, including the transparent or semitransparent fluid flows in high-speed motion, usually outperforms the capability of traditional monitoring, a strain gauge system consists by dual-plane sensors, EMF flow meter (OPTI-FLUX 4000, from KROHNE), and pressure sensors was used as an effective supplementary for flow measurement (Liang, Ye, Wang, & Brauwer, 2012). It was mounted on the backside of sheet-typed ring by 2-D array with vertical and span wise axes, to inspect the shock-wave vibration and stress distribution  of waterjet energy beam concurrently. All the obtained signals were finally sent to a 64-bit data acquisition board connected. As the grinding conditions were exampled by Table 2, a total of 100 measurements were conducted for each workpiece specimen. After establishing a multidimensional matrix array for data arrangement, the measured planar distribution of turbulence characteristics, and the specific coordinate positions can be processed accurately or trans-dimensional highlighted by vector diagrams (Pellerin, Leclaire, & Reggio, 2015). Some typical examples of characteristic demonstration were shown in Figures 2 and 5. Based on the actual measured turbulence characteristics by LDA and strain gauge systems, a general picture describing the turbulence characteristics can be established (Nadaoka, Nihei, & Yagi, 1999).  Table 11 presented their mean values at different time points, and provided a standard benchmark for the prediction evaluations with typical flow models. It was clearly seen that as the AWJ grinding continues, the turbulence characteristics in the contact zone presented a more scattered and decentralized tendency due to the localized stepping on the machined surface.
With the help of those actual-measured values, orthogonal experiment design was employed, thereafter the number of experimental tests can be greatly decreased, offering an optimal factor combination for actual usage. The following steps of orthogonal experiment were conducted in sequence: Determining the experiment objective (turbulence characteristics) and establishing the target function (fuzzy prediction function); Selecting the influence indexes (AWJ working parameters); Deciding the representative value levels for objective indexes; Establishing a suitable orthogonal array; Analyzing the CFD simulations or machining experiments; And finally determining an optimum combination for influence factors and value levels (Hattori, Hotta, & Houra, 2014). For example, Table 12 presents the value levels of AWJ parameters, Table 13 proposes the partitioned universes for turbulence characteristics in Reynolds stress model. Using an orthogonal array of experimental scheme for the waterjet pressure, traverse speed, flow velocity, abrasive flow rate, and attack angle yielded 40 test combinations. Table 14 gives some representative orthogonal rules for the prediction verification, Table 15 presents the coding identification and the comparison of characteristic results.
To give an accurate assessment on the effectiveness of orthogonal experiment, F-ratio-tests were employed to classify the total variations of data values into two parts: the variation caused by the change of working parameter, and the variation caused by the experimental error (Hiziroglu, 2013). The F-ratio value demonstrates that the targeted influence factor exerts a greater impact on the fuzzy prediction of turbulence characteristics than other factor does. In addition, the contribution rates of working parameters were quantified to identify their relative importance (Sarkheyli, Zain, & Sharif, 2015).
The F-ratio criteria was shown as follows: w p denotes the number of working parameters, there were m j levels for the j th parameter, each level was tested for c times; n denotes the total number of orthogonal tests; the obtained value for each verification case came from the measure points was marked as x 1 , x 2 , . . . , x n , K 1 , K 2 , . . . , K mj responded to the level sum of the j th parameter with respect to the testing times (Shabgard et al., 2013): Thus, the square sum of variance: SS j = 1 c m j i=1 K i 2 − (1/m j c)K 2 ; the degree of variance freedom: df j = m j -1; The total square sum: SS T = n i=1 x i 2 − (1/m j c)K 2 ; the total degree of variance freedom: df T = n-1; The square sum of error: SS e = SS T − w p j=1 SS j ; the degree error of variance freedom: df e = df T − w p j=1 df j ; The mean square value S j = SS j df j ; the mean square × value of error MS e = SS e df e ; F -ratio The importance analysis of targeted factors was shown by Table 16. It can be observed that the effects of P w , W a and T s in k-model were obviously. P w and F a presented a greater impact than other factors did in k--R t flow model. Simultaneously, the influence of W a in Spalart-Allmaras model cannot be neglected. Further investigation shows that P w and T s exerted powerful influences    when LES model was used, while W a or F a become crucial in Reynolds stress model (Greifzu, Kratzsch, Forgber, Lindner, & Schwarze, 2016). Deviations between the predicted and actual measured characteristics can be learned in Figures 16-20. In these figures, the horizontal axis denotes the predicted turbulence characteristic, while the vertical axis gives the actual measured one, the ratios between them by using different flow models were demonstrated by colored dots, and emphasized by black circles.
In order to compare the prediction performances influenced by different flow models, several benchmark coefficients were investigated to assess their actual working capabilities in an identical experimental condition: the response time (T r : The time consumption for characteristic prediction and was measured by automated time-counting instruments), the occupied memory storage (M s : The occupied memory capacity used for fuzzy  1  1  1  1  1  1  2  2  4  1  3  2  1  1  1  2  2  5  1  2  2  2  3  1  1  2  3  3  3  2  5  5  5  prediction and was calibrated from the real-time performance indicator) (Pakhomov, & Terekhov, 2015;Zhang, Law, & Zhao, 2015). Furthermore, some mathematical coefficients have been proposed to assess the prediction accuracy in practice (Wang, Tang, & Wu, 2016) Here EV and PV denote the predicted and actual measured values, the meanings of other variables can be learned form the literature (Liang, Liu, Ye, & Wang, 2013); As the statistical metric results were offered at the bottom of Figures 16-20, which supports the influence investigation of typical flow models as respected. Figure 16 presents the ratios between the predicted and actual measured RMS velocities. In this case it can be learned that k--R t model outperforms others in the prediction accuracy, for its ratio dots scattered close around the straight line of slope 45°where the accuracy of characteristic prediction was calibrated. Correspondingly, the ratio dots denoting k-and Spalart-Allmaras models kept above the straight line of slope 45°, demonstrating their remarkable influences upon the fuzzy prediction of RMS velocity in a negative way. On the other hand, the ratio dots of LES and Reynolds stress model situated below the straight line, as a result of over-prediction in RMS velocity. Since the energy dissipation and pressure gradient reached an optimum state when k--R t model was employed, its widely application has been acknowledged in RMS velocity prediction already.
Turbulence intensity ( Figure 17) kept a high reliability when Spalart-Allmaras model was employed. The computation results of k-model showed an obviously increasing tendency, and deviating from the benchmark line thereafter. k--R t model generated an obviously fluctuation in the actual fuzzy prediction, as the black dots demonstrated. LES model can be freely used to show the distribution of turbulence intensity on workpiece surface, simultaneously Reynolds stress model gave an important impact on the fuzzy prediction of turbulence intensity.
The fuzzy prediction of turbulence kinetic energy ( Figure 18) kept a closely correlation with Spalart-Allmaras model in the proposed experimental conditions. It can be easily observed that the ratio dots denoted by Spalart-Allmaras model kept a high agreement with the predetermined straight line, which verifying its accuracy and reliability in kinetic energy prediction even when the external noise signal or work-hardening capacities were introduced in. In an actual grinding condition, this flow model presents a reliable criterion to assess the predicted outcome of kinetic energy. In a contrast, the fuzzy prediction of k--R t or Reynolds stress model cannot eliminate the external interferences resulted from turbulence disturbance or error prediction, therefore they cannot ensure a clearer result demonstration. The ratio dots of LES or k-model showed that they can make a robust prediction only when the turbulence strain was considered about in waterjet grinding.
In the prediction of turbulence entropy (Figure 19), k--R t model enjoys its unique priority in evaluation accuracy, especially when the measurement deviation was focused on. From repeatedly experiments and detailed data analysis it was learned that k-and LES models presented a less accurate prediction of entropy when a radical flow vortex happened. The prediction comparison between the Reynolds stress and Spalart-Allmaras models shows that, although some improvements were realized to minimize the prediction tolerance and truncation error, the deviation error between the predicted and actual-measured results still unsolved unless the frequently characteristic sampling and appropriate computation tolerances were accompanied with. k--R t model obtains a precise illustration of turbulence entropy when the stress tensor or energy dissipation on the boundary layer were emphasized, that explains its widely The k-model used for the prediction of Reynolds shear stress (Figure 20) returns an accurate result, as the amendment quantity of turbulence stress loading maintains stable during the whole waterjet grinding process. The prediction comparison between k--R t and LES models demonstrates that the predicted results remained in a relatively unstable state, consequently the computation precision was impacted negatively (Yao, Yang, & Wang, 2016;Zanaganeh, Mousavi, & Etemad Shahidi, 2009). In some interferential conditions such as turbulence shock and radical waterjet impact, Spalart-Allmaras and Reynolds stress models can provide a relatively stable information feedback or standard signal benchmark for turbulence monitoring.
For the purpose of assessing the performances of fuzzy prediction by using different flow models in a visual way, Figures 21-25 described the overall analysis of them in radar diagrams for clearer illustration and better understanding, with seven diagram axles stand for the selected benchmark coefficients mentioned before: T r , M s , MAE, MAPE, NMAE, NRMSE and Pearson_r. It can be found that when V RMS was fuzzy predicted, k-model ensures a good performance in T r , MAPE, NMAE, and NRMSE, since it was more suitable than others to demonstrate the distribution of V RMS throughout the targeted section of waterjet stream, and especially appropriate to give an objective assessment of dynamic fluid properties in the contact zone; k--R t model returns an excellent characteristic prediction of T i when the NRMSE or Pear-son_r were based upon, since it realizes an accurate  description of data relationship in the turbulence property sets, and achieves a remarkable development in T r , M s or MAE when compared with that of other models; Spalart-Allmaras model provides a more accurate prediction in T c than others, characterized by the criteria of MAPE, NMAE, and NRMSE. This model gives a straightforward prediction of flow property distribution in the statistical domain, including some relevant turbulence characteristics such as T e , V RMS , or R s . Reynolds stress model owns a good working capability in the fuzzy prediction of T e , which provided a more robust prediction result in the view of NMAE, NRMSE and Pearson_r, than other alternatives in the waterjet machining process, thus a helpful technical reference in T r and M s can be  provided, too (Rybalko, Loth, & Lankford, 2012;); Finally, the excellent performances of LES model outperformed that of other flow models in R s , M s , MAPE and NMAE, demonstrating the fact that this model provides a more accurate and comprehensive prediction result of R s in actual machining practices. Similar working capabilities or performance evaluations can be learned in the cases of T i , T c , and T e (Wu,Chau,& Fan,Figure 21. The performance comparisons of fuzzy prediction for V RMS by using different flow models.    2010; Wu, Chau, & Li, 2008). After summarizing their respective prediction results, it is noteworthy to point out that the prediction accuracy of flow models based upon the precondition that the complicated interaction between the working conditions and turbulence characteristics was an empirical relationship, therefore, their applications and superiorities can be verified definitely.

Conclusions
After the fuzzy predictions and performance comparisons, a set of integrated assessments and helpful suggestions can be made on the specific applications of flow models for turbulence characteristic predictions: k--R t model should be widely used for the preliminary characteristic analysis of waterjet flow; k-model was usually used to predict the planar distribution of Reynolds shear stress and turbulence intensity in the contact area; Spalart-Allmaras model can be conveniently employed for turbulence characteristic assessments by monitoring the pressure gradient or fluid velocity; LES model shows its superiorities in the turbulence intensity description when the impulsive loading was exerted; Reynolds stress model can be applied to obtain a reliable prediction result of turbulence characteristic when faced with the highly energy dissipation, stress tensor, or vortex momentum. All these assessments and suggestions provided valuable references and practical instructions to predict the AWJ grinding performance. Furthermore, (1) As other studies have not realized the AWJ flow description in the contact zone, a series of mathematical characteristics describing the multiphase waterjet flow were proposed quantitatively; (2) Different from other traditional methods concluding the turbulence characteristics from a macro-scale dimensional simulation, a new fuzzy prediction system has been proposed to predict them, and distinguish the influences of typical flow models on the turbulence characteristic predictions; (3) As the ordinary studies simply focused on the predictions of AWJ grinding performance without any further considerations of the important influences caused by flow models, this research tried to make comparisons between the fuzzy prediction results of AWJ turbulence characteristics, by using orthogonal experiment with the CFD simulations and actual characteristic measurements; (4) The final conclusions were obtained based upon the discussions of influence mechanism and theoretical superiorities of typical flow models.
Since an universal description concentrating on the influence effectiveness of flow models absents and becomes a major difficult problem in the current waterjet grinding researches, developing a more robust prediction system and making a comprehensive influence investigation for the turbulence characteristic in AWJ grinding conditions, will be heart of our future efforts, therefore more machining comparisons and theoretical analysis between the multi-phase flow models were suggested, in the interests of obtaining more reliable predictions of AWJ machining effectiveness. Additionally, it was suggested that the correlations between the waterjet characteristics and the machining performances be studied, to provide the technical references for AWJ grinding monitoring in actual conditions.

Disclosure statement
No potential conflict of interest was reported by the authors.