INTRODUCTION
Nanofluid firstly introduced by Dr. Choi in Argonne National Laboratory describes a fluid in which nanometersize particles are suspended (Choi et al., 2001). Nanoparticles can be highly useful for various industrial applications due to their unique properties such as large surface area to volume ratio and lower kinematic energy (Yu et al., 2007; Keblinski et al., 2005). Graphene as one of the most efficient nanoparticles has a highquality electron transport at room temperature. Various theoretical studies have been performed so far to determine the thermal conductivity of graphene particles (Saito et al., 2007). For example, in an experimental study carried out by Balandin et al., thermal conductivity of a singlelayer graphene was calculated whose value was found to be 5300 W/mK (Balandin et al., 2008). Due to the aforementioned characteristics, nanofluids are extensively utilized in electronic and cooling industries. In a research reported by Theres et al., thermal conductivities and heat transfer properties of several different nanofluids containing Hydrogen exfoliated graphene (HEG), dispersed deionized (DI) water, and ethylene glycol (EG) were systematically investigated (Baby and Ramaprabhu, 2011). They observed that increment in Nusselt number was much more than that in thermal conductivity for the considered nanofluids. Keblinski et al. described four possible mechanisms to justify the unique increment of thermal conductivity in nanofluids (Paul et al., 2012): a) molecularlevel layering of the liquid at liquid/particle interface; b) Brownian motion of nanoparticles; c) the nature of heat transport in nanoparticles; and d) the effect of nanoparticle clustering. Brownian motion of nanoparticles is usually too slow to conduct thermal energy through the nanofluid directly. However, this motion can play an important indirect role in inducing convection and clustering of nanoparticles which results in heat transfer enhancement (Azimi and Ommi, 2013). Nowadays, identification of squeezing flow characteristics has attracted much attention among scholars due to its importance in improving the performance of hydraulic machine elements and numerous applications in food industry, chemical engineering, polymer processing, compression and injection molding (Ghori et al., 2007; Ran et al., 2009; Burbidge and Servais, 2004; Grimm, 1976; Siddiqui et al., 2008). Since unexpected behavior of fluid flows can lead to undesirable repercussions like mechanical noise in engineering systems, analysis of the behaviour of the flows under different conditions is indispensable (Jianu and Rosen, 2017). It is obvious that almost all of scientific phenomena such as heat and mass transfer processes occur nonlinearly. In fluid mechanics problems, inherent nonlinear behavior of the systems should be considered for reaching accurate and precise results. Solving nonlinear equations are often very complicated and cannot be done by ordinary methods. Sheikholeslami et al. analyzed an unsteady nanofluid flow between two parallel plates by using (Adomain Decomposition Method) ADM (Sheikholeslami et al., 2013). The results showed that when both the plates are moveable, Nusselt number declines with an increase in squeeze number and increases with increments in nanoparticle volume fraction and Eckert number. Sheikholeslami and Ganji analytically examined the hydrothermal behavior of nanofluid flows taking into account magnetic field and Brownian motion by using Differential Transformation Method (Sheikholeslami and Ganji, 2015). The results indicated that skin friction coefficient rises with increases in squeeze and Hartmann numbers and decreases with an increment in nanofluid volume fraction. In addition, they concluded that Nusselt number increases with increasing the nanoparticle volume fraction and Hartmann number while it decreases with growing the squeeze number. Dogonchi et al. used an analytical solution known as DuanRach approach to study the flow and heat transfer of nanofluid within an unsteady squeezing flow considering the thermal radiation effect (Dogonchi et al., 2016). According to their results, temperature profile and Nusselt number increase with rising the radiation parameter. Azimi and Mozaffari considered an unsteady two dimensional Graphene Oxide water based nanofluid to examine the heat energy transferred between two moving parallel plates utilizing an intelligent blackbox identifier (Azimi and Mozaffari, 2015). Based on their results, it is inferred that the developed (hybrid genetic mutable smart bee algorithmfuzzy inference system) HGMSBAFIS blackbox identifier can be applied as an authentic tool for analyzing the flow due to its accuracy and robustness. Azimi et al. considered an unsteady nanofluid flow between two moving parallel plates and analytically solved the nonlinear governing equations using Galerkin Optimal Homotopy Asymptotic Method (GOHAM) (Azimi et al., 2014). The results obtained from their numerical analysis via fourth grade order Runge–Kutta had satisfactory agreement with the results yielded from GOHAM. Sheikholeslami and Ganji used Homotopy Perturbation Method to analyze the thermal behavior of a nanofluid squeezing flow (Sheikholeslami and Ganji, 2013). With respect to their results, it can be extracted that Nusselt number has direct relationship with nanoparticle volume fraction, squeeze and Eckert numbers when the two plates are separate and has reverse relationship with squeeze number when the two plates are squeezed. RahimiGorji et al. utilized Galerkin method to solve nonlinear differential equations in an unsteady squeezing nanofluid flow and heat transfer problem based on Brinkman model considering variable magnetic field (RahimiGorji et al., 2016). According to their results, as nanofluid volume fraction increases, Nusselt number rises and skin friction coefficient decreases. Moreover, enhancement in nanofluid volume fraction increases the Hartman number resulting in an increase in velocity and temperature profiles. Recently, Sheikholeslami and Ganji conducted a numerical analysis on CuOwater nanofluid to describe the heat transfer variations occurring between two coaxial tubes in presence of magnetic field (Sheikholeslami and Ganji, 2017). Their results showed that temperature gradient increases with an increment in Hartmann number and decreases with an increase in melting parameter. Furthermore, velocity reduces with increasing the Lorentz forces and melting parameter while it rises with an increment in Reynolds number. In another study carried out by Sheikholeslami et al., an analytical solution is developed for a nanofluid flow to investigate the heat transfer between two pipes in existence of magnetic field using AGM method (Sheikholeslami et al., 2017). The results indicated that temperature gradient enhances with an increase in Hartman and Eckert numbers declines with an enhancement in Reynolds number.
In this study, the temperature behavior of an unsteady flow containing water as base fluid and graphene nanoparticles between two moving plates in a twodimensional configuration is examined for the first time using AkbariGanji method. Viscous dissipation is considered in the governing equations and reliabl functions are presented for temperature and velocity distributions using AGM. Finally, the impact of viscous dissipation are elaborately explained for several nanoparticles and Eckert numbers.
GOVERNING EQUATIONS
As previously mentioned, in this study, an unsteady nanofluid flow and heat transfer in a twodimensional configuration between two infinite parallel plates is considered and investigated. Mixture of liquid water (base fluid) and Graphene Oxide (nanoparticle) is taken as the nanofluid. Schematic of the problem considered in the current analysis is indicated in Figure 1. The viscous dissipation and heat generation parameters caused by friction and shear stress are considered in the present study. These parameters play important roles particularly in large values of Eckert number. Eckert number describes the relationship between the flow kinetic energy and enthalpy. Following assumptions are used for the analysis of the nanofluid flow:

The flow is assumed to be incompressible.

Nochemical reaction occurs in the system.

The effect of radiative heat transfer is negligible.

Solid nanoparticles and base fluid (water) are in thermal equilibrium.

Noslip condition exists between the solid nanoparticles and water.
Timedependent position of the plates can be calculated by the following correlation (Sheikholeslami et al., 2013):
\[z = \pm {l(1  \alpha t)}^{\frac{1}{2}} = \pm h(t)\]
Governing equations can be expressed as follows for the problem (Sheikholeslami et al., 2013):
\[\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0\] 
(1) 
\[\rho_{\text{nf}}\left( \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} \right) =  \frac{\partial p}{\partial x} + \mu_{\text{nf}}\left( \frac{\partial^{2}u}{\partial x^{2}} + \frac{\partial^{2}u}{\partial y^{2}} \right)\] 
(2) 
\[\rho_{\text{nf}}\left( \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} \right) =  \frac{\partial p}{\partial y} + \mu_{\text{nf}}\left( \frac{\partial^{2}v}{\partial x^{2}} + \frac{\partial^{2}v}{\partial y^{2}} \right)\] 
(3) 
\[\left( \frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} \right) = \frac{k_{\text{nf}}}{{(\rho C_{p})}_{\text{nf}}}\left( \frac{\partial^{2}T}{\partial x^{2}} + \frac{\partial^{2}T}{\partial y^{2}} \right) + \frac{\mu_{\text{nf}}}{{(\rho C_{p})}_{\text{nf}}}\left\lbrack 4{(\frac{\partial u}{\partial x})}^{2} + {(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y})}^{2} \right\rbrack\] 
(4) 
where \(u\) and \(v\) are velocities in x and y directions, respectively. For the nanofluid,\(\ \rho_{\text{nf}}\), \(\mu_{\text{nf}}\), \({(\rho C_{p})}_{\text{nf}}\) and \(k_{\text{nf}}\) are effective density, effective dynamic viscosity, effective heat capacity and effective thermal conduction coefficient, respectively, and can be defined as follows (Sheikholeslami et al., 2013):
\[\left\{ \begin{matrix} \rho_{\text{nf}} = \rho_{f}(1  \varphi) + \rho_{s}\varphi_{s} \\ {(\rho C_{p})}_{\text{nf}} = \left( \rho C_{p} \right)_{f}\left( 1  \varphi \right) + {(\rho C_{p})}_{s} \\ \mu_{\text{nf}} = \frac{\mu_{f}}{\left( 1  \varphi \right)^{2.5}} \\ \frac{k_{\text{nf}}}{k_{f}} = \frac{k_{s} + 2k_{f}  2\varphi\left( k_{f}  k_{s} \right)}{k_{s} + 2k_{f} + 2\varphi\left( k_{f}  k_{s} \right)} \\ \vartheta_{\text{nf}} = \frac{\mu_{f}}{\rho_{\text{nf}}} \\ \end{matrix} \right.\ \] 
(5) 
It should be noted that in the above equations, thermal conductivity coefficient and viscosity are calculated by using MaxwellGarnetts (MG) and Brinkman models, respectively.
Boundary conditions are expressed as follows (Sheikholeslami et al., 2013):
\[\frac{y = h\left( t \right) \rightarrow v = v_{w} = \frac{\text{dh}}{\text{dt}},\ \ T = T_{H}}{y = 0 \rightarrow v = \frac{\partial u}{\partial y} = \frac{\partial T}{\partial y} = 0}\] 
(6) 
In order to nondimensionalize the governing equations, following parameters are defined (Sheikholeslami et al., 2013):
\[\frac{\eta = \frac{y}{\left\lbrack l{(1  \alpha t)}^{\frac{1}{2}} \right\rbrack},\ u = \frac{\text{αx}}{2\left( 1  \alpha t \right)}f^{'}\left( \eta \right),\ \ \ v = \frac{\text{αl}}{\left\lbrack \left( 1  \alpha t \right)^{\frac{1}{2}} \right\rbrack}f(\eta)}{\theta = \frac{T}{T_{H}},\ A = \left( 1  \varphi \right) + \varphi\frac{\rho_{s}}{\rho_{f}},\ \ B = \left( 1  \varphi \right) + \varphi\frac{({\rho C_{p})}_{s}}{{(\rho C_{p})}_{f}},\ C = \frac{k_{\text{nf}}}{k_{f}}}\] 
(7) 
where \(\alpha(t)\) and \(v_{w} = \frac{d\alpha(t)}{\text{dt}}\) are the distance of each plate from x coordinate and velocity of the plates, respectively. When \(\alpha > 0\), with increasing the time, the plates get away from each other and when \(\alpha < 0\), the plates approach each other which is known as squeezing flow.
By eliminating the pressure gradient term from the governing equations and applying the aforementioned nondimensionalized parameters presented in Eq. 7, conservation laws are expressed as follows (Sheikholeslami et al., 2013):
\[f^{(IV)}  SA\left( 1  \varphi \right)^{2.5}\left( \eta f^{'''} + 3f^{''} + f^{'}f^{''}  ff^{''} \right) = 0\] 
(8) 
\[\theta^{''} + PrS\left( \frac{B}{C} \right)\left( f\theta^{'}  \eta\theta^{'} \right) + \frac{\Pr\text{EC}}{C\left( 1  \varphi \right)^{2.5}}({f^{''}}^{2} + 4\delta^{2}{f^{'}}^{2}) = 0\] 
(9) 
Considering the nondimensionalized parameters, following boundary conditions will be achieved (Sheikholeslami et al., 2013):
\[f\left( 0 \right) = 0,\ f^{''}\left( 0 \right) = 0,\ f\left( 1 \right) = 1,\ f^{'}\left( 1 \right) = 0,\ \ \ \theta^{'}\left( 0 \right) = 0,\ \ \ \theta\left( 1 \right) = 1\] 
(10) 
where \(S\), \(\Pr\) and \(\text{Ec}\) are moving parameter, prandtl number and Eckert number, respectively. The abovementioned parameters can be defined as follows (Sheikholeslami et al., 2013):
\[S = \frac{\alpha l^{2}}{2v_{f}},\ \ \ Pr = \frac{\mu_{f}{(\rho C_{p})}_{f}}{\rho_{f}k_{f}},\ \ \ Ec = \frac{\rho_{f}}{{(\rho C_{p})}_{f}}{(\frac{\text{αx}}{2(1  \alpha t)})}^{2},\ \delta = \frac{1}{x}\] 
(11) 
If \(S > 0\), the plates get away from each other and if \(S < 0\) the plates move toward each other.
AKBARIGANJI METHOD (AGM)
According to this method, general form of an equation (Eq. 12) with its boundary conditions (Eq. 13) are defined as follows (Sheikholeslami et al., 2017):
\[P_{k}:f\left( u,u^{'},u^{''},\ldots,u^{m} \right) = 0\ ;\ \ \ \ \ u = u(x)\] 
(12) 
\[\left\{ \begin{matrix} u\left( 0 \right) = u_{0}, u^{'}\left( 0 \right) = u_{1},\ldots,u^{m  1}\left( 0 \right) = u_{m  1} \\ u\left( l \right) = u_{l_{0}}, u^{'}\left( l \right) = u_{l_{1}},\ldots,u^{m  1}\left( l \right) = u_{l_{m  1}} \\ \end{matrix} \right.\ \] 
(13) 
It is assumed that following equation is the solution of this problem (Sheikholeslami et al., 2017):
\[u\left( x \right) = \sum_{i = 0}^{n}{a_{i}x^{i} = a_{0} + a_{1}x^{1} + a_{2}x^{2} + \ldots + a_{n}x^{n}}\] 
(14) 
The higher the values of n, the more the accuracy of the solution. By inserting Eq. (14) into Eq. (12), the residual is obtained. Regarding the boundary conditions and the values of residual at the boundaries, constant parameters of Eq. (14) are calculated.
RESULTS AND DISCUSSION
As previously mentioned, the nonlinear governing equations for the flow between two moving parallel plates are solved using AkbariGanji method. The derivations and calculations are conducted in Matlab environment. Table 1 lists the properties of the considered fluid (water) and nanoparticles. It is worth mentioning that the current numerical analysis is undertaken using Boundary Value Problem (BVP) function.
Table 1. Thermophysical properties of water and the nanoparticles
Fluid

\[\mathbf{\rho} \] (kg/m3)

\[\mathbf{C_p}\] (j/kg.K)

\[\mathbf{K}\] (W/m.K)

Pure water

997.1

4179

0.613

Graphene Oxide

1800

717

5000

Copper (Cu)

8933

385

401

Silver (Ag)

10500

235

429

Alumina (Al_{2}O_{3})

3970

765

40

Titanium Oxide (TiO_{2})

4250

686.2

8.9538


Validation
Regarding the defined boundary conditions, following approximate functions are considered. It is notable that these functions are obtained and selected based on try and error approach and their performance in the problem.
\[f = {a_{0} + a_{1}\eta + a_{2}\eta^{2} + a_{3}\eta}^{3} + a_{4}\eta^{4} + a_{5}\eta^{5} + a_{6}\eta^{6} + a_{7}\eta^{7} + a_{8}\eta^{8}\] 
(15) 
\[\theta = {b_{0} + b_{1}\eta + b_{2}\eta^{2} + b_{3}\eta}^{3} + b_{4}\eta^{4} + b_{5}\eta^{5} + b_{6}\eta^{6} + b_{7}\eta^{7} + b_{8}\eta^{8}\] 
(16) 
With respect to these functions and by enforcing the methods, unknown parameters will be calculated. In this study, nondimensionalized forms of volume fraction (\(\delta\)), plate moving parameter, Prandtl number and Eckert number are assumed to be 0.05, 0.1, 1, 6.2 and 0.1, respectively. Following functions will be achieved by inserting Eqs. (15) and (16) into the nondimensionalized governing equations and doing some algebraic calculations:
\[f = {1.4621\eta  0.5\eta}^{3} + 0.1673\eta^{4}  0.1697\eta^{5} + 0.0594\eta^{6}  0.0132\eta^{7}  0.0057\eta^{8}\] 
(17) 
\[\theta = {1.3864  0.0247\eta^{2}  0.4149\eta}^{4} + 0.345\eta^{5}  0.8449\eta^{6} + 1.0355\eta^{7}  0.4825\eta^{8}\] 
(18) 
Figure 2 illustrates the results of the numerical and analytical solutions conducted for solving function f. As can be seen in Figure 2, there is a good agreement between the results of the current analytical and numerical solutions conducted on momentum equation and therefore AkbariGanji method can be an efficient and desirable method for solving momentum equation.
Results of the numerical and analytical solutions for function f are listed in Table 2. As can be observed in Table 2, the values of errors between the results of the numerical analysis and analytical solution are very small. According to the table, maximum and average values of error are found to be about 1.3 % and 0. 71 %, respectively.
Table 2. Comparison between the errors of the numerical analysis and analytical solution for function f
𝞰

Error of AGM

0

0

0.1

0.004963

0.2

0.0091

0.3

0.011872

0.4

0.013004

0.5

0.012462

0.6

0.010449

0.7

0.007405

0.8

0.004016

0.9

0.001191

1

4.29E09


Results of the numerical and the analytical analyses performed for function \(\theta\) are delineated in Figure 3. Regarding Figure 3, it is observed that results of the analytical solution and the numerical analysis are very close to each other for energy equation. In this regard, in addition to the momentum equation, this method gives us appropriate and acceptable results for energy equation. By comparing Figures 2 and 3, it can be implied that errors related to the solution of energy equation is a bit more than those related to the calculation of momentum equation. In order to improve the precision of the results for energy equation, the temperature function can be approximated by a polynomial function with higher degrees. Since in the present investigation the values of errors are small, the first order of the function is used.
Errors related to the numerical and analytical solutions for calculating the energy equation are indicated in Figure 3. As it can be seen in this figure, there is an acceptable agreement between the numerical and analytical results and this method can be an efficient approach for solving energy equation. According to Figure 4, maximum and average values of error are found to be 3% and 1.59%, respectively.
Effect of viscous dissipation on temperature distribution
In this subsection, the effect of viscous dissipation on temperature distribution is presented. In energy differential equation (Eq. 9), the last term on the lefthand side is related to the viscous dissipation which can be calculated after obtaining the function f and twice integrating\(\ \frac{\Pr\text{EC}}{C\left( 1  \varphi \right)^{2.5}}({f^{''}}^{2} + 4\delta^{2}{f^{'}}^{2})\). Afterwards, the influence of the viscous dissipation on the temperature distribution is examined. In this investigation, the ratio of the result achieved by integrating to the total temperature obtained from Eq. (9) is introduced as viscous dissipation ratio. As can be implied from this definition, Eckert number is an effective parameter in this ratio. This number generally represents the ratio of kinetic energy to the enthalpy difference within the boundary layer and its importance is mainly in highvelocity flows and in considerable values of viscous dissipation. Figure 5 shows the effect of position on the viscous dissipation ratio. As previously mentioned, nondimensionalized forms of volume fraction (\(\delta\)), plate moving parameter, Prandtl number and Eckert number are assumed to be 0.05, 0.1, 1, 6.2 and 0.1, respectively. With respect to Figure 5, viscous dissipation ratio has larger values near the plates. According to this figure, with an increment in the value of position, the magnitude of this ratio decreases pointing out that the impact of viscosity is highly effective near the plates.
Figure 6 shows the influence of Eckert number on viscous dissipation ratio. Regarding Figure 6, an increase in the value of Eckert number results in an increment in the value of viscous dissipation ratio and amount of heat energy caused by shear stress. Variation of temperature with position is depicted in Figure 7 for several Eckert numbers. As can be observed in this figure, with increasing the value of Eckert number, temperature varies further and the amount of heat transfer increases.
Figure 8 indicates the mean values of viscous dissipation ratio for several types of nanoparticles within the flow domain. Based on Figure 8, it can easily be concluded that the average value of viscous dissipation ratio of titanium oxide and graphene is considerably more than that of the other considered nanoparticles. It is worth noting that the viscous dissipation ratio of the nanoparticles, for instance silver and titanium oxide, are very close to each other as shown in Figure 9. In other words, changing the type of nanoparticles has negligible effect on the viscous dissipation ratio. In order to illustrate the insignificant differences in Figure 9 better and with more precision, zoomin view of results are presented in Figure 10.
Effect of squeeze number on Velocity and temperature profiles
Figure 11 presents the effect of squeeze number on the velocity. It should be stated that the squeeze number (S) describes the movement of the plates (S > 0 corresponds to the plates moving apart, while S < 0 corresponds to the plates moving together. The positive and negative squeeze numbers have different impacts on the velocity profile. Velocity increases with an enhancement in the absolute value of squeeze number when η < 0.5 and decreases when η > 0.5. It is worth noting that when η= 0.5, velocity is constant and moving parameter has no effect on the velocity. Figure 12 delineates the impact of squeeze number on the temperature. When the two plates move toward each other, thermal boundary layer thickness grows with increasing the absolute magnitude of the squeeze number. This increment in thermal boundary layer thickness leads to a reduction in Nusselt number. Moreover, an opposite behavior is observed when two plates get away from each other. Consequently, for cooling applications, higher values of squeeze number is desirable and for heating applications, the lower values of squeeze number is favorable.
CONCLUSION
In the present study, a twodimensional unsteady nanofluid flow containing water as the base fluid and graphene nanoparticles between two moving parallel plates was studied. Governing partial differential equations were initially turned into ordinary differential equations utilizing similarity solution. Energy and momentum equations were analytically solved by AkbariGanji method and reliable mathematical functions were obtained for temperature and velocity distributions. Moreover, in order to check the accuracy of the method, the governing equations were also solved by numerical solutions. Satisfactory agreement was observed between the analytical and numerical results. In this regard, it can be concluded that the method can be efficient for solving nonlinear ordinary differential equations. The values of errors associated with the momentum and energy equations were in orders of 0.001 and 0.01, respectively. Furthermore, viscous dissipation ratio was calculated for several nanoparticle types and Eckert numbers. Results indicated that with increasing the Eckert number, the effect of viscous dissipation ratio on generated heat caused by shear stress would be further. Ultimately, the impact of nanoparticle type on the viscous dissipation was investigated. According to the results, changing the type of nanoparticle had negligible effect was on the viscous dissipation and when titanium oxide nanoparticles was used, maximum value of viscous dissipation ratio was achieved.
NOMENCLATURE
\[u\] 
Velocity in x direction 
\[v\] 
Velocity in y direction 
\[p\] 
Pressure 
\[\rho_{\text{nf}}\] 
Effective density of nanofluid 
\[{{(\text{ρC}}_{P})}_{\text{nf}}\] 
Effective heat capacity of nanofluid 
\[k_{\text{nf}}\] 
Effective thermal conductivity of nanofluid 
\[\mu_{\text{nf}}\] 
Viscosity of nanofluid 
\[\Pr\] 
Prandtl number 
\[S\] 
Moving parameter 
Greek symbols

\[\alpha\] 
Constant rotational velocity 
\[\varphi\] 
Dimensionless concentration 
\[\mu\] 
Dynamic viscosity 
\[\vartheta\] 
Kinematic viscosity 
\[\theta\] 
Dimensionless temperature 
\[\rho\] 
Fluid density 
Subscripts 
\[f\] 
Base fluid 
\[\text{nf}\] 
Nanofluid 
\[p\] 
Nano particle 