Nomenclature
| $C$ | fluid concentration |
| ${C}_{\infty }$ | concentration at infinity |
| ${C}_{g}$ | gel layer concentration |
| ${\tilde{D}}_{\gamma }$ | generalized diffusion coefficient |
| ${J}_{y}$ | mass transfer flux |
| $k$ | permeability parameter |
| $L$ | characteristic length |
| ${n}_{1},{n}_{2}$ | constants |
| $\mathop{R{\rm{e}}}\limits^{\sim }$ | generalized Reynolds number |
| $\mathop{Sc}\limits^{\sim }$ | generalized Schmidt number |
| $t$ | time |
| $u,$ $v$ | velocity components along the $x$ and $y$ directions, respectively |
| ${u}_{\infty }$ | free-flow velocity |
| ${v}_{w}$ | penetration velocity |
| $x,$ $y$ | coordinates along the plate and normal to it, respectively |
| Greek symbols | |
| $\alpha $ | fractional spatial derivative of velocity |
| $\gamma $ | fractional spatial derivative of concentration |
| ${\tilde{\mu }}_{\alpha }$ | generalized dynamic viscosity |
| $\rho $ | fluid density |
| ${\tau }_{yx}$ | shear stress |
1. Introduction
Filtration is widely utilized in numerous industrial processes. For instance, in liquid food processing, the fluid flowing through a penetrable surface is always accompanied by suction/injection. Owing to its wide range of applications, surface penetration has become a topic of considerable research interest. Fakhar et al [1] studied rotational flow on a penetrable plate filled with a second-order fluid using Lie symmetry analysis. Moatimid et al [2] investigated the electrodynamic Kelvin–Helmholtz instability when viscous potential fluid flowing between the interface of two porous layers, and the impact of a horizontal electric field and suction/injection were considered. Maity [3] addressed the flow and heat transfer problems of liquid films with suction/injection through a penetrable stretching sheet. Aziz et al [4] examined the unsteady flow of a third-grade fluid flowing in a porous medium on a permeable sheet. The Soret–Dufour effect on the heat and mass transfer of non-Newtonian fluids through a perpendicular penetrable surface was discussed by Idowu et al [5]. The flow driven by a rotating disk under the impact of uniform injection/suction in a porous media was analyzed by Attia [6], who also examined the heat transfer. Considering the influence of suction/injection, Sheremet et al [7] examined the mixed convection of a nanofluid within a square porous cavity. Sheikholeslami [8] addressed the natural convection of a nanofluid within a permeable enclosure, and the existence of Lorentz forces was considered. Furthermore, the effect of an external magnetic source on the heat transfer of a Fe3O4-water nanofluid filled in a permeable enclosure was analyzed by Sheikholeslami and Seyednezhad [9]. Sheikholeslami et al [10] investigated the forced convection of nanofluids filled in a porous enclosure by using the Lattice Boltzmann method, and some thermal circular obstacles were considered. Ellahi et al [11] addressed magnetohydrodynamic flow of non-Newtonian fluids in saturating porous spaces. Hassan et al [12] analyzed the micropolar flow over a permeable channel with mass injection. Considering the chemical reaction and thermal absorption/generation effect, Eid et al [13] analyzed the boundary layer of a Carreau-type nanofluid flowing along a nonlinear stretching sheet immersed in a porous media. Arasteh et al [14] performed a numerical simulation on nanofluids flowing in a three-layer porous medium and investigated the heat transfer characteristics.
In practical industrial applications, it is more common to observe non-uniform suction/injection, and many studies have been conducted on this subject. Saikrishnan et al [15] examined the free convection and non-uniform mass transfer over a cylinder or a sphere, considering the non-uniform slot blowing/suction effect. Roy et al [16] analyzed the boundary layer flowing on a rotating sphere under the impacts of heterogeneous slot injection/suction. Furthermore, they addressed the boundary layer flow past an inclined cylinder, and the suction/injection on the surface was observed to be variable [17]. Kumari et al [18] focused on Newtonian flow along a cone vertically inserted in a saturated porous medium, and it was assumed that the suction/injection followed a power-law variation with the distance. Ganapathirao et al [19] addressed the boundary layer flow adjacent to a wedge, and non-uniform suction/injection was considered. Ravindran et al [20] studied the steady mixed convection of viscous and incompressible fluids flowing along a vertical cone with non-uniform suction/injection effect. Li et al [21] analyzed the flow and heat transfer of a power-law liquid food over a non-uniform penetrable channel. Bai et al [22] addressed the flow, heat, and mass transfer problem of Maxwell fluid past a bidirectional stretching sheet. The non-uniform suction/injection velocity of fluid flowing through the porous permeable plate was written as a specific function ${v}_{w}(x,t)=k{t}^{{n}_{1}}{x}^{{n}_{2}}{y}^{{n}_{3}},$ where $k,$ ${n}_{1},$ ${n}_{2}$ and ${n}_{3}$ are constants. We have adopted a function similar to that in [22] to describe the non-uniform permeation in this study. Furthermore, Mondal et al [23] studied the rheological model of non-Newtonian fluid migration in the ultrafiltration process of a juice. In the concentration boundary layer, they considered steady-state mass balance, and the boundary condition ${v}_{w}{C}_{g}+{\left.D\tfrac{\partial C}{\partial y}\right|}_{y=0}=0$ was established, which has been referenced in our research.
In recent years, fractional constitutive models have attracted researchers' attention owing to their advantages in depicting ‘time and distance memory'. Sheikh et al [24] studied the natural convection of a magnetohydrodynamic second-order fluid filled in porous media using a time-fractional derivative model. Carrera et al [25] derived a fractional Maxwell model to describe the viscous behavior of non-Newtonian fluid. Sun et al [26] proposed a new model of spatial fractional-order constitutive equation that showed a relation between shear stress and fractional-order velocity gradient, and it could suitably depict complex non-Newtonian fluid. In [27], the Riemann–Liouville fractional derivative was used by Pan et al to relate the stress relaxation of viscoelastic fluids with distance. Li et al [28] addressed the flow of viscoelastic fluid through a non-uniform permeable surface by using the fractional constitutive model. Furthermore, scholars applied fractional constitutive models to liquid food. According to the experimental results of [29], the fractional-order model can be successfully used to describe the viscoelastic creep behavior of a Hami melon. Simpson et al [30] proposed that the diffusion in food materials was abnormal and that fractional calculus could be used to describe the mass transfer in food processing. Inspired by above studies, the Riemann–Liouville fractional spatial derivative is used in this study to describe the stress relaxation and mass transfer of viscoelastic liquid food. In addition, the concentration boundary conditions are modified according to the fractional calculus model, and a new boundary condition ${v}_{w}(x,t){C}_{g}+\mathop{{D}_{\gamma }}\limits^{\sim }\tfrac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }}(x,0,t)=0$ of mass transfer is established.
Fractional differential equations have been widely adopted in the fields of anomalous diffusion and hydrodynamics. However, deriving the analytic solutions of fractional differential equations has been difficult. Consequently, several numerical algorithms have been proposed and tested. For example, the classical Galerkin finite-element method was adopted by Mehdi et al [31] to solve the spatial fractional diffusion-wave equation defined by Riemann–Liouville fractional derivative. Liu et al [32] solved the nonlinear fourth-order diffusion equation numerically with a Galerkin method which used second-order time approximation scheme, and the stability of the method was proved. Yang et al [33] used the Legendre spectrum method in the space direction and L1 finite difference scheme in the time direction to solve the first Stokes problem of a generalized fractional second-order fluid. Ding et al [34] used higher-order numerical approximation scheme to solve the spatial fractional convection equation, and analyzed the stability of the algorithm. Yang et al [35] constructed a novel method with higher accuracy and efficiency by modifying the coefficients of the upwind scheme, and the effective conditions applicable to the scheme has been analyzed. Feng et al [36] solved a series of two-dimensional sub-diffusion and diffusion-wave equations with multiple time-fractional differential terms on convex domains via the finite difference/element method. The numerical solution of the generalized Maxwell fluid flow and heat transfer were obtained by Liu et al [37] using temporal fractional constitutive models with L1 and L2 schemes. Hao et al [38] proposed the quasi-compact finite difference scheme, and they solved the space fractional equation with delay numerically and analyzed the convergence of the scheme. Sayevand et al [39] solved the fractional Navier–Stokes equations via a new non-standard finite difference algorithm. Liu et al [40] dealt with the space-time-fractional advection-dispersion equation using the finite difference method, and they proved the stability and convergence of the algorithm. Zhang et al [41] solved the governing equations of Maxwell nanofluids with the time and space-fractional-order conditions using the finite difference method. The rotating electro-osmotic flow of fractional Maxwell fluids was addressed by Wang et al [42] utilizing finite difference method. We adopted the finite difference method herein and tested the accuracy and mesh dependence of the algorithm.
Motivated by the abovementioned studies, the flow and mass transfer process of viscoelastic liquid food during filtration through a permeable plate is investigated. In section 2 , the governing boundary layer equations depicting fluid flow and mass transfer are established based on fractional constitutive models. In section 3 , the Grünwald–Letnikov approximation discretization schemes of the Riemann–Liouville fractional derivative are adopted, and the finite difference algorithm is used to solve the problem. In section 4 , the accuracy of the algorithm and the influence of step size on the results are verified. In section 5 , the influences of time, fractional constitutive models, generalized Reynolds number, generalized Schmidt number, and the permeability parameter on the velocity and concentration fields are discussed.
2. Mathematical formulation
In this study, the flow of a viscoelastic fluid (e.g. low-in-fat oil-in-water food emulsions [43], corn starch-milk-sugar paste [44]) is investigated through a permeable plate, as shown in figure 1. During food processing, the liquid food flows through a permeable surface, long enough to ensure that the fluid flow is fully developed. The fluid has a uniform mainstream velocity ${u}_{\infty }.$ On the surface of the plate, during the food filtration process, the pores will be blocked by the settled impurities in the fluid, or the liquid food will spontaneously form a gel layer deposition [23]. Thus, the fluid will permeate at an uneven rate ${v}_{w}(x,t).$ It is assumed that the fluid concentration is uniform at the inlet, and the fluid density is constant during the process of penetration. In order to ensure the effect of suction/injection does not interference the flow field in the boundary layer, the penetration velocity ${v}_{w}(x,t)$ is chosen to be relative small. Furthermore, chemical reactions are negligible, and the concentration of the gel layer is below the solubility limit; thus, the continuity equation still applies.
Figure 1. Schematic of viscoelastic incompressible liquid food flowing through a permeable plate. |
The governing equations are given as follows.
Continuity equation:
$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial x}+\displaystyle \frac{\partial v}{\partial y}=0.\end{eqnarray}$
Momentum and concentration equations:
$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial t}+u\displaystyle \frac{\partial u}{\partial x}+v\displaystyle \frac{\partial u}{\partial y}=\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial {\tau }_{yx}}{\partial y},\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{\partial C}{\partial t}+u\displaystyle \frac{\partial C}{\partial x}+v\displaystyle \frac{\partial C}{\partial y}=\displaystyle \frac{\partial {J}_{y}}{\partial y}.\end{eqnarray}$
It is well known that viscoelastic liquid food is usually a mixture, which is composed of substances of different types such as water, oil, polymers and other long chain molecules. The motion and dynamic characteristics of these mixtures often exhibit long-term memory and spatial nonlocal properties [45]. However, Newton's constitutive model is characterized by integral-order velocity gradient, which only reveals the local effect of the system. The fractional derivative model, which has the properties of spatial nonlocality and long-term memory, has recently been selected as a useful tool of describing viscoelastic fluids. For instance, the flow of viscoelastic fluid in a channel was characterized successfully by Qi et al [46] using the fractional Maxwell model. The fractional constitutive model has been used by Yang et al [47] to study the start-up flow of a viscoelastic fluid in a pipe. Ding et al [48] characterized hydraulic fracturing in viscoelastic formation adopting the fractional model. In this paper, we investigate the boundary layer of the viscoelastic fluid, of which the shear stress at a certain point is not only dependent on the velocity gradient of this point but also the velocity gradient of points nearby or even far away. As in [27], we adopt a fractional-order relation between the shear stress and the velocity gradient:
$\begin{eqnarray}{\tau }_{yx}={\tilde{\mu }}_{\alpha }\displaystyle \frac{{\partial }^{\alpha }u}{\partial {y}^{\alpha }},\end{eqnarray}$
where ${\tilde{\mu }}_{\alpha }$ denotes the generalized dynamic viscosity, and $\alpha $ $(0\lt \alpha \lt 1)$ is the space-fractional derivative. The Riemann–Liouville fractional derivative is used [49]: $\begin{eqnarray}\displaystyle \frac{{\partial }^{\alpha }u}{\partial {y}^{\alpha }}=\displaystyle \frac{1}{{\rm{\Gamma }}(1-\alpha )}\displaystyle \frac{\partial }{\partial y}\displaystyle {\int }_{0}^{y}{(y-\xi )}^{-\alpha }u(x,\xi ,t){\rm{d}}\xi ,\end{eqnarray}$
where ${\rm{\Gamma }}(\cdot )$ denotes the Gamma function.Furthermore, abnormal diffusion is common in food and biomaterials processing. Welti-Chanes et al [50] summarized that in many cases of food processing, non-Fickian behaviors existed in convection phenomenon caused by buoyancy or tissue expansion, stirring effect, and active/passive membrane transport. In the process of osmotic dehydration of many vegetables and fruits, the internal mass transfer depended on the microstructure of biological tissues, and abnormal diffusion appeared [51]. Simpson et al [52] investigated the diffusion mechanism of apple osmotic dehydration in moderate electric field, and it showed that the sucrose diffusion in apple slices did not follow the conventional diffusion model except in 17 V cm−1 electric field. Ramírez et al [53] demonstrated that fractional calculus could relatively well characterize the abnormal diffusion in apple slices. Watanabe et al [54] confirmed that the water transport in starchy foods did not follow the conventional diffusion model. Furthermore, salting is also one of the common processes in food processing, where abnormal diffusion can be observed. By simulating the diffusion of salt and water in salmon halogenation process, Núñez et al [55] demonstrated that the anomalous model based on fractional calculus could better describe the diffusion phenomenon. For a viscoelastic food liquid, the mass transfer flux at a certain point can also be affected by the concentration gradient of points nearby or even those far away. Similar to [30], the fractional-order relation between mass transfer flux and concentration gradient is employed:
$\begin{eqnarray}{J}_{y}={\tilde{D}}_{\gamma }\displaystyle \frac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }},\end{eqnarray}$
where ${\tilde{D}}_{\gamma }$ denotes the generalized diffusion coefficient, and $\gamma $ $(0\lt \gamma \lt 1)$ is the fractional spatial derivative. The Riemann–Liouville fractional derivative is also adopted as in [49]: $\begin{eqnarray}\displaystyle \frac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }}=\displaystyle \frac{1}{{\rm{\Gamma }}(1-\gamma )}\displaystyle \frac{\partial }{\partial y}\displaystyle {\int }_{0}^{y}{(y-\xi )}^{-\gamma }C(x,\xi ,t){\rm{d}}\xi .\end{eqnarray}$
By using the fractional-order relations of (4 ) and (6 ), the governing equations (2 ) and (3 ) become
$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial t}+u\displaystyle \frac{\partial u}{\partial x}+v\displaystyle \frac{\partial u}{\partial y}=\displaystyle \frac{{\tilde{\mu }}_{\alpha }}{\rho }\displaystyle \frac{\partial }{\partial y}\left(\displaystyle \frac{{\partial }^{\alpha }u}{\partial {y}^{\alpha }}\right),\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{\partial C}{\partial t}+u\displaystyle \frac{\partial C}{\partial x}+v\displaystyle \frac{\partial C}{\partial y}={\tilde{D}}_{\gamma }\displaystyle \frac{\partial }{\partial y}\left(\displaystyle \frac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }}\right).\end{eqnarray}$
The initial and boundary conditions are as follows:12 ). ${C}_{g}$ is the concentration of the gel layer.
$\begin{eqnarray}\begin{array}{l}u(x,y,0)={u}_{\infty },\,v(x,y,0)=0,\,C(x,y,0)={C}_{\infty },\\ x\geqslant 0,\,y\gt 0,\end{array}\end{eqnarray}$
$\begin{eqnarray}u(x,0,t)=0,\,v(x,0,t)={v}_{w}(x,t),\,x\gt 0,\,t\gt 0,\end{eqnarray}$
$\begin{eqnarray}{v}_{w}(x,t){C}_{g}+\mathop{{D}_{\gamma }}\limits^{\sim }\displaystyle \frac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }}(x,0,t)=0,\,x\gt 0,\,t\gt 0,\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u(x,+\infty ,t)={u}_{\infty },\,v(x,+\infty ,t)=0,\\ C(x,+\infty ,t)={C}_{\infty },\,x\geqslant 0,\,t\gt 0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u(0,y,t)={u}_{\infty },\,v(0,y,t)=0,\,C(0,y,t)={C}_{\infty },\\ y\gt 0,\,t\geqslant 0,\end{array}\end{eqnarray}$
where ${v}_{w}(x,t)$ is the rate of penetration and can be written as a function of distance and time, i.e. ${v}_{w}(x,t)=k{x}^{{n}_{1}}{t}^{{n}_{2}}$ as in [22], where $k,$ ${n}_{1}$ and ${n}_{2}$ are constants. This function is employed to describe a two-dimensional problem and reveal the fact that in the process of infiltration, impurities in the fluid precipitate or accumulate, resulting in an uneven inner wall of the filter tube and a non-uniform fluid penetration. It is noted that $k$ is the permeability parameter, when $k\lt 0$ corresponds to suction, and $k\gt 0$ represents injection. Furthermore, in [23], Mondal et al considered the steady-state mass balance during ultrafiltration on the plate surface, and the boundary condition ${v}_{w}{C}_{g}+{\left.D\tfrac{\partial C}{\partial y}\right|}_{y=0}=0$ was established. The first item on the left side of the equation denotes the mass transfer flux caused by penetration, while the second item is the mass transfer flux caused by molecular diffusion. When abnormal diffusion happens in food and biomaterials processing, the conventional Fick's law $J=D\tfrac{\partial C}{\partial y}$ cannot better describe the mass transfer flux caused by molecular diffusion. Thus, the modified Fick's law ${J}_{y}={\tilde{D}}_{\gamma }\tfrac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }}$ is employed and a fractional-order boundary condition makes the problem real. Our study improved the concentration boundary conditions according to the modified Fick's law and obtained equation (The following dimensionless variables are adopted:
$\begin{eqnarray}\begin{array}{l}x* \,=\,\displaystyle \frac{x}{L},\,y* \,=\,\displaystyle \frac{y}{L},\,u* \,=\,\displaystyle \frac{u}{{u}_{\infty }},\,v* \,=\,\displaystyle \frac{v}{{u}_{\infty }},\,{v}_{w}^{\ast }=\displaystyle \frac{{v}_{w}}{{u}_{\infty }},\\ C* \,=\,\displaystyle \frac{C}{{C}_{\infty }},\,{C}_{g}^{\ast }=\displaystyle \frac{{C}_{g}}{{C}_{\infty }},\,t* \,=\,\displaystyle \frac{t{u}_{\infty }}{L},\,k* \,=\,\displaystyle \frac{{L}^{{n}_{1}+{n}_{2}}}{{u}_{\infty }^{{n}_{2}+1}}k,\end{array}\end{eqnarray}$
where characteristic length is denoted by $L.$ The dimensionless governing equations become $\begin{eqnarray}\displaystyle \frac{\partial u* }{\partial x* }+\displaystyle \frac{\partial v* }{\partial y* }=0,\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{\partial u* }{\partial t* }+u* \displaystyle \frac{\partial u* }{\partial x* }+v* \displaystyle \frac{\partial u* }{\partial y* }=\displaystyle \frac{1}{\tilde{R{\rm{e}}}}\displaystyle \frac{\partial }{\partial y* }\left(\displaystyle \frac{{\partial }^{\alpha }u* }{\partial {y}^{\ast \alpha }}\right),\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{\partial C* }{\partial t* }+u* \displaystyle \frac{\partial C* }{\partial x* }+v* \displaystyle \frac{\partial C* }{\partial y* }=\displaystyle \frac{1}{\tilde{R{\rm{e}}}}\displaystyle \frac{1}{\tilde{Sc}}\displaystyle \frac{\partial }{\partial y* }\left(\displaystyle \frac{{\partial }^{\gamma }C* }{\partial {y}^{\ast \gamma }}\right),\end{eqnarray}$
where $\mathop{R{\rm{e}}}\limits^{\sim }=\tfrac{\rho {u}_{\infty }{L}^{\alpha }}{{\tilde{\mu }}_{\alpha }}$ is the generalized Reynolds number, and $\mathop{Sc}\limits^{\sim }=\tfrac{{\tilde{\mu }}_{\alpha }}{\rho {L}^{\alpha -\gamma }{\tilde{D}}_{\gamma }}$ is the generalized Schmidt number.The corresponding non-dimensional initial and boundary conditions are
$\begin{eqnarray}\begin{array}{l}u* (x* ,\,y* ,\,0)=1,v* (x* ,\,y* ,\,0)=0,\\ C* (x* ,\,y* ,\,0)=1,x* \,\geqslant \,0,y* \,\gt \,0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u* (x* ,\,0,\,t* )=0,v* (x* ,\,0,\,t* )={v}_{w}^{\ast }(x* ,\,t* ),\\ x* \,\gt \,0,t* \,\gt \,0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\mathop{R{\rm{e}}}\limits^{\sim }\mathop{Sc}\limits^{\sim }{C}_{g}^{\ast }{v}_{w}^{\ast }(x* ,\,t* )+\displaystyle \frac{{\partial }^{\gamma }C* }{\partial {y}^{\ast \gamma }}(x* ,\,0,t* )=0,\\ x* \,\gt \,0,t* \,\gt \,0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u* (x* ,+\infty ,\,t* )=1,v* (x* ,\,+\infty ,t* )=0,\\ C* (x* ,+\infty ,t* )=1,x* \,\geqslant \,0,t* \gt 0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u* (0,y* ,\,t* )=1,v* (0,\,y* ,\,t* )=0,\\ C* (0,y* ,\,t* )=1,y* \,\gt \,0,t* \,\geqslant \,0,\end{array}\end{eqnarray}$
where $\begin{eqnarray}{v}_{w}^{\ast }(x* ,\,t* )=k* {x}^{* {n}_{1}}{t}^{* {n}_{2}}.\end{eqnarray}$
The governing equations (16 )–(18 ) and boundary conditions (19 )–(23 ) are to be solved.
3. Implicit numerical method
Equations (16 )–(18 ) are calculated by decoupling and linearization as follows.
First, the general form of equations (17 ), (18 ) is re-written (for simplicity we delete mark* hereafter) as2 , $f(x,y,t)$ and $g(x,y,t)$ are both equal to $0.$ $a(x,y,t)$ corresponds to the horizontal velocity $u(x,y,t),$ which is greater than or equal to 0, and $b(x,y,t)$ denotes the vertical velocity $v(x,y,t).$ According to equation (16 ), $a(x,y,t)$ and $b(x,y,t)$ satisfy the following formula, i.e.
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial u}{\partial t}+a(x,y,t)\displaystyle \frac{\partial u}{\partial x}+b(x,y,t)\displaystyle \frac{\partial u}{\partial y}={F}_{\alpha }\displaystyle \frac{{\partial }^{\alpha +1}u}{\partial {y}^{\alpha +1}}\\ \,\,\,+\,f(x,y,t),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial C}{\partial t}+a(x,y,t)\displaystyle \frac{\partial C}{\partial x}+b(x,y,t)\displaystyle \frac{\partial C}{\partial y}={L}_{\gamma }\displaystyle \frac{{\partial }^{\gamma +1}C}{\partial {y}^{\gamma +1}}\\ \,\,\,\,+\,g(x,y,t),\end{array}\end{eqnarray}$
where ${F}_{\alpha }=\tfrac{1}{\mathop{R{\rm{e}}}\limits^{\sim }},$ ${L}_{\gamma }=\tfrac{1}{\mathop{R{\rm{e}}}\limits^{\sim }}\tfrac{1}{\mathop{Sc}\limits^{\sim }}.$ In the physical problem presented in section $\begin{eqnarray}\displaystyle \frac{\partial a(x,y,t)}{\partial x}+\displaystyle \frac{\partial b(x,y,t)}{\partial y}=0.\end{eqnarray}$
Define ${t}_{k}=k\tau ,$ $k=0,1,\,\ldots ,\,n,$ ${x}_{i}={\rm{i}}{h}_{x},$ $i=0,1,\,\ldots ,{M}_{x},$ ${y}_{j}=j{h}_{y},$ $j=0,1,\,\ldots ,\,{M}_{y},$ where $\tau =T/n$ is time step; ${h}_{x}$ is the space step in $x$ direction; and ${h}_{y}$ is the space step in $y$ direction. Assume $u(x,y,t)\in {C}^{2},$ $v(x,y,t)\in {C}^{2},$ $C(x,y,t)\in {C}^{2}.$
The Grünwald–Letnikov approximation discretization schemes of the Riemann–Liouville derivatives ${\partial }^{\alpha +1}u/\partial {y}^{\alpha +1}$ and ${\partial }^{\gamma +1}C/\partial {y}^{\gamma +1}$ are as follows [56]:
$\begin{eqnarray}\begin{array}{l}\frac{{\partial }^{\alpha +1}}{\partial {y}^{\alpha +1}}u({x}_{i},{y}_{j},{t}_{k+1})=\frac{1}{{h}_{y}^{\alpha +1}}\displaystyle \sum _{m=0}^{j+1}{\omega }_{m}^{\alpha +1}u\left({x}_{i},{y}_{j+1-m},\right.\\ \,\,\,\,\,\left.{t}_{k+1}\right)+O({h}_{y}),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\frac{{\partial }^{\gamma +1}}{\partial {y}^{\gamma +1}}C({x}_{i},{y}_{j},{t}_{k+1})=\frac{1}{{h}_{y}^{\gamma +1}}\displaystyle \sum _{q=0}^{j+1}{\eta }_{q}^{\gamma +1}C\left({x}_{i},{y}_{j+1-q},\right.\\ \,\,\,\,\,\left.{t}_{k+1}\right)+O({h}_{y}).\end{array}\end{eqnarray}$
Taking as
$\begin{eqnarray*}\begin{array}{l}{\omega }_{0}^{\alpha +1}=1,\,{\omega }_{m}^{\alpha +1}=\left(1-\displaystyle \frac{\alpha +2}{m}\right){\omega }_{m-1}^{\alpha +1},\,m=1,2,\ldots ,\\ {\eta }_{0}^{\gamma +1}=1,\,{\eta }_{q}^{\gamma +1}=\left(1-\displaystyle \frac{\gamma +2}{q}\right){\eta }_{q-1}^{\gamma +1},\,q=1,2,\ldots .\end{array}\end{eqnarray*}$
Using the backward difference method, $\partial u/\partial x$ and $\partial C/\partial x$ can be approximated as
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial u({x}_{i},{y}_{j},{t}_{k+1})}{\partial x}=\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i-1},{y}_{j},{t}_{k+1})}{{h}_{x}}\\ \,+\,O({h}_{x}),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial C({x}_{i},{y}_{j},{t}_{k+1})}{\partial x}=\displaystyle \frac{C({x}_{i},{y}_{j},{t}_{k+1})-C({x}_{i-1},{y}_{j},{t}_{k+1})}{{h}_{x}}\\ \,+\,O({h}_{x}).\end{array}\end{eqnarray}$
Similarly, $\partial u/\partial y$ and $\partial C/\partial y$ are approximated as
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial u({x}_{i},{y}_{j},{t}_{k+1})}{\partial y}=\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i},{y}_{j-1},{t}_{k+1})}{{h}_{y}}\\ \,+\,O({h}_{y}),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial C({x}_{i},{y}_{j},{t}_{k+1})}{\partial y}=\displaystyle \frac{C({x}_{i},{y}_{j},{t}_{k+1})-C({x}_{i},{y}_{j-1},{t}_{k+1})}{{h}_{y}}\\ \,+\,O({h}_{y}).\end{array}\end{eqnarray}$
For the first-order temporal derivative, the approximate discrete scheme is obtained by using the backward difference scheme:
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial u({x}_{i},{y}_{j},{t}_{k+1})}{\partial t}=\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i},{y}_{j},{t}_{k})}{\tau }\\ \,+\,O(\tau ),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial C({x}_{i},{y}_{j},{t}_{k+1})}{\partial t}=\displaystyle \frac{C({x}_{i},{y}_{j},{t}_{k+1})-C({x}_{i},{y}_{j},{t}_{k})}{\tau }\\ \,+\,O(\tau ).\end{array}\end{eqnarray}$
Hence, equations (25 ) and (26 ) become25 ) utilizing the implicit difference scheme is as follows:
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i},{y}_{j},{t}_{k})}{\tau }=-a({x}_{i},{y}_{j},{t}_{k})\\ \,\times \,\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i-1},{y}_{j},{t}_{k+1})}{{h}_{x}}-b({x}_{i},{y}_{j},{t}_{k})\\ \,\times \,\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i},{y}_{j-1},{t}_{k+1})}{{h}_{y}}\\ \,+\,\displaystyle \frac{{F}_{\alpha }}{{h}_{y}^{\alpha +1}}\displaystyle \sum _{m=0}^{j+1}{\omega }_{m}^{\alpha +1}u({x}_{i},{y}_{j+1-m},{t}_{k+1})\\ \,+\,f({x}_{i},{y}_{j},{t}_{k+1})+{R}_{i,j,k+1},\end{array}\end{eqnarray}$
and $\begin{eqnarray}\begin{array}{l}\displaystyle \frac{C({x}_{i},{y}_{j},{t}_{k+1})-C({x}_{i},{y}_{j},{t}_{k})}{\tau }=-a({x}_{i},{y}_{j},{t}_{k+1})\\ \,\times \,\displaystyle \frac{C({x}_{i},{y}_{j},{t}_{k+1})-C({x}_{i-1},{y}_{j},{t}_{k+1})}{{h}_{x}}\\ \,-\,b({x}_{i},{y}_{j},{t}_{k+1})\displaystyle \frac{C({x}_{i},{y}_{j},{t}_{k+1})-C({x}_{i},{y}_{j-1},{t}_{k+1})}{{h}_{y}}\\ \,+\,\displaystyle \frac{{L}_{\gamma }}{{h}_{y}^{\gamma +1}}\displaystyle \sum _{q=0}^{j+1}{\eta }_{q}^{\gamma +1}C({x}_{i},{y}_{j+1-q},{t}_{k+1})\\ \,+\,g({x}_{i},{y}_{j},{t}_{k+1})+{Z}_{i,j,k+1},\end{array}\end{eqnarray}$
where $\begin{eqnarray}\left|{R}_{i,j,k+1}\right|\leqslant A(\tau +{h}_{x}+{h}_{y}),\end{eqnarray}$
$\begin{eqnarray}\left|{Z}_{i,j,k+1}\right|\leqslant B(\tau +{h}_{x}+{h}_{y}).\end{eqnarray}$
$A$ and $B$ are constants. Note that ${u}_{i,j}^{k+1}$ is a numerical approximation of $u({x}_{i},{y}_{j},{t}_{k+1}).$ Then, the approximation of equation ( $\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{u}_{i,j}^{k+1}-{u}_{i,j}^{k}}{\tau }=-{a}_{i,j}^{k}\displaystyle \frac{{u}_{i,j}^{k+1}-{u}_{i-1,j}^{k+1}}{{h}_{x}}-{b}_{i,j}^{k}\displaystyle \frac{{u}_{i,j}^{k+1}-{u}_{i,j-1}^{k+1}}{{h}_{y}}\\ \,\,\,\,\,+\,\displaystyle \frac{{F}_{\alpha }}{{h}_{y}^{\alpha +1}}\displaystyle \sum _{m=0}^{j+1}{\omega }_{m}^{\alpha +1}{u}_{i,j+1-m}^{k+1}+{f}_{i,j}^{k+1},\end{array}\end{eqnarray}$
or $\begin{eqnarray}\begin{array}{l}{\mu }_{i,j}^{k}{u}_{i,j}^{k+1}-{r}_{1}{a}_{i,j}^{k}{u}_{i-1,j}^{k+1}-{r}_{2}{b}_{i,j}^{k}{u}_{i,j-1}^{k+1}-{r}_{3}\displaystyle \sum _{m=0,m\ne 1}^{j+1}{\omega }_{m}^{\alpha +1}\\ \,\,\,\,\times \,{u}_{i,j+1-m}^{k+1}={u}_{i,j}^{k}+\tau {f}_{i,j}^{k+1},\end{array}\end{eqnarray}$
where ${r}_{1}=\tfrac{\tau }{{h}_{x}},$ ${r}_{2}=\tfrac{\tau }{{h}_{y}},$ ${r}_{3}=\tfrac{\tau {F}_{\alpha }}{{h}_{y}^{\alpha +1}},$ ${\mu }_{i,j}^{k}=1+{r}_{1}{a}_{i,j}^{k}\,+{r}_{2}{b}_{i,j}^{k}-{r}_{3}{\omega }_{1}^{\alpha +1},$ ${a}_{i,j}^{k}$ is the numerical approximation of $u({x}_{i},{y}_{j},{t}_{k}),$ ${b}_{i,j}^{k}$ is the numerical approximation of $v({x}_{i},{y}_{j},{t}_{k}).$Similar as above, ${C}_{i,j}^{k+1}$ is the numerical approximation of $C({x}_{i},{y}_{j},{t}_{k+1}).$ The approximation of equation (26 ) is as follows:
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{C}_{i,j}^{k+1}-{C}_{i,j}^{k}}{\tau }=-{a}_{i,j}^{k+1}\displaystyle \frac{{C}_{i,j}^{k+1}-{C}_{i-1,j}^{k+1}}{{h}_{x}}-{b}_{i,j}^{k+1}\\ \,\,\,\,\,\,\,\times \,\displaystyle \frac{{C}_{i,j}^{k+1}-{C}_{i,j-1}^{k+1}}{{h}_{y}}+\displaystyle \frac{{L}_{\gamma }}{{h}_{y}^{\gamma +1}}\displaystyle \sum _{q=0}^{j+1}{\eta }_{q}^{\gamma +1}{C}_{i,j+1-q}^{k+1}+{g}_{i,j}^{k+1},\end{array}\end{eqnarray}$
or $\begin{eqnarray}\begin{array}{l}{\lambda }_{i,j}^{k+1}{C}_{i,j}^{k+1}-{r}_{1}{a}_{i,j}^{k+1}{C}_{i-1,j}^{k+1}-{r}_{2}{b}_{i,j}^{k+1}{C}_{i,j-1}^{k+1}-{r}_{4}\displaystyle \sum _{q=0,q\ne 1}^{j+1}{\eta }_{q}^{\gamma +1}\\ \,\,\,\,\,\,\times \,{C}_{i,j+1-q}^{k+1}={C}_{i,j}^{k}+\tau {g}_{i,j}^{k+1},\end{array}\end{eqnarray}$
where ${r}_{1}=\tfrac{\tau }{{h}_{x}},$ ${r}_{2}=\tfrac{\tau }{{h}_{y}},$ ${r}_{4}=\tfrac{\tau {L}_{\gamma }}{{h}_{y}^{\gamma +1}},$ ${\lambda }_{i,j}^{k+1}=1+{r}_{1}{a}_{i,j}^{k+1}\,+{r}_{2}{b}_{i,j}^{k+1}-{r}_{4}{\eta }_{1}^{\gamma +1},$ ${a}_{i,j}^{k+1}$ is the numerical approximation of $u({x}_{i},{y}_{j},{t}_{k+1}),$ ${b}_{i,j}^{k+1}$ is the numerical approximation of $v({x}_{i},{y}_{j},{t}_{k+1}).$We now consider the boundary condition (21 ). The Grünwald–Letnikov approximation discretization schemes for the Riemann–Liouville derivative ${\partial }^{\gamma }C/\partial {y}^{\gamma }$ are also adopted, i.e.
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{\partial }^{\gamma }}{\partial {y}^{\gamma }}C({x}_{i},{y}_{j},{t}_{k+1})=\displaystyle \frac{1}{{h}_{y}^{\gamma }}\displaystyle \sum _{p=0}^{j+1}{\chi }_{p}^{\gamma }C({x}_{i},{y}_{j+1-p},{t}_{k+1})\\ \,+\,O({h}_{y}).\end{array}\end{eqnarray}$
Taking as ${\chi }_{0}^{\gamma }=1,$ ${\chi }_{p}^{\gamma }=\left(1-\tfrac{\gamma +1}{p}\right){\chi }_{p-1}^{\gamma },$ $p=1,2,\,\ldots .$
Boundary condition (21 ) becomes
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{1}{{h}_{y}^{\gamma }}\displaystyle \sum _{p=0}^{j+1}{\chi }_{p}^{\gamma }C({x}_{i},{y}_{j+1-p},{t}_{k+1})=-\tilde{{Re}}\tilde{Sc}{C}_{g}{v}_{w}({x}_{i},{t}_{k+1})\\ \,+\,{M}_{i,j,k+1},j=1,\end{array}\end{eqnarray}$
where $\begin{eqnarray}\left|{M}_{i,j,k+1}\right|\leqslant N(\tau +{h}_{x}+{h}_{y}).\end{eqnarray}$
$N$ is a constant.
Adopting ${C}_{i,j}^{k+1}$ as the approximation of $C({x}_{i},{y}_{j},{t}_{k+1}),$ the implicit difference approximation of boundary condition (21 ) is as follows:43 ) and (48 ), the following equation is obtained:
$\begin{eqnarray}{d}_{1}\displaystyle \sum _{p=0}^{j+1}{\chi }_{p}^{\gamma }{C}_{i,j+1-p}^{k+1}=-\tilde{B}{v}_{w\,i}^{k+1},\,j=1,\end{eqnarray}$
or $\begin{eqnarray}{C}_{i,0}^{k+1}=\left(-\displaystyle \frac{\tilde{B}}{{d}_{1}}{v}_{w\,i}^{k+1}-{\chi }_{0}^{\gamma }{C}_{i,1}^{k+1}\right)/{\chi }_{1}^{\gamma },\end{eqnarray}$
where $\tilde{B}=\mathop{{Re}}\limits^{\sim }\mathop{Sc}\limits^{\sim }{C}_{g},$ ${d}_{1}=\tfrac{1}{{h}_{y}^{\gamma }}.$ Combining equations ( $\begin{eqnarray}\begin{array}{l}{\lambda }_{i,j}^{k+1}{C}_{i,j}^{k+1}-{r}_{1}{a}_{i,j}^{k+1}{C}_{i-1,j}^{k+1}-{r}_{2}{b}_{i,j}^{k+1}{C}_{i,j-1}^{k+1}-{r}_{4}\\ \,\times \displaystyle \sum _{q=0,q\ne 1}^{j}{\eta }_{q}^{\gamma +1}{C}_{i,j+1-q}^{k+1}-{r}_{4}{\eta }_{j+1}^{\gamma +1}{C}_{i,0}^{k+1}={C}_{i,j}^{k}+\tau {g}_{i,j}^{k+1},\end{array}\end{eqnarray}$
i.e. $\begin{eqnarray}\begin{array}{l}{\lambda }_{i,j}^{k+1}{C}_{i,j}^{k+1}-{r}_{1}{a}_{i,j}^{k+1}{C}_{i-1,j}^{k+1}-{r}_{2}{b}_{i,j}^{k+1}{C}_{i,j-1}^{k+1}-{r}_{4}\displaystyle \sum _{q=0,q\ne 1}^{j}{\eta }_{q}^{\gamma +1}\\ \,\times \,{C}_{i,j+1-q}^{k+1}+{r}_{4}{\eta }_{j+1}^{\gamma +1}\displaystyle \frac{{\chi }_{0}^{\gamma }}{{\chi }_{1}^{\gamma }}{C}_{i,1}^{k+1}={C}_{i,j}^{k}+\tau {g}_{i,j}^{k+1}-{\theta }_{i,j}^{k+1},\end{array}\end{eqnarray}$
where ${\theta }_{i,j}^{k+1}={r}_{4}{\eta }_{j+1}^{\gamma +1}\tfrac{\tilde{B}}{{d}_{1}{\chi }_{1}^{\gamma }}{v}_{w\,i}^{k+1}.$Note that $v({x}_{i},{y}_{j},{t}_{k+1})$ is obtained from equation (27 ):
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{u({x}_{i},{y}_{j},{t}_{k+1})-u({x}_{i-1},{y}_{j},{t}_{k+1})}{{h}_{x}}+O({h}_{x})\\ \,=-\displaystyle \frac{v({x}_{i},{y}_{j},{t}_{k+1})-v({x}_{i},{y}_{j-1},{t}_{k+1})}{{h}_{y}}+O({h}_{y}).\end{array}\end{eqnarray}$
Similar as above, ${v}_{i,j}^{k+1}$ is the numerical approximation of $v({x}_{i},{y}_{j},{t}_{k+1}).$ The approximation of equation (27 ) is as follows:
$\begin{eqnarray}{v}_{i,j}^{k+1}=\displaystyle \frac{{h}_{y}}{{h}_{x}}({u}_{i-1,j}^{k+1}-{u}_{i,j}^{k+1})+{v}_{i,j-1}^{k+1}.\end{eqnarray}$
4. Algorithm analysis
4.1. Mesh independence
By calculating the governing equations (16 )–(18 ) and boundary conditions (19 )–(23 ), the effect of step size on the results is compared, and the appropriate step size for subsequent work is selected.
The computed region is treated as a rectangle. In figures 2–4, the spatial step ${h}_{x},$ ${h}_{y}$ and time step $\tau $ varies to demonstrate their influences on the velocity $u$ and concentration $C.$ The results prove that the proposed finite-difference algorithm is stable when the grid sizes are limited above certain numbers.
Figure 2. Comparison of $u$ and $C$ with different space step size ${h}_{x}$ with $x=0.8,$ $t=1,$ $\alpha =0.8,$ $\gamma =0.7,$ $\tilde{R{\rm{e}}}=4,$ $\tilde{Sc}=4,$ ${C}_{g}=6.18,$ $k=-5\times {10}^{-3},$ ${n}_{1}=2,$ ${n}_{2}=3,$ ${h}_{y}=1/50,$ $\tau =1/50.$ |
Figure 3. Comparison of $u$ and $C$ with different space step size ${h}_{y}$ with $t=1,$ $\alpha =0.8,$ $\gamma =0.7,$ $\mathop{R{\rm{e}}}\limits^{\sim }=4,$ $\mathop{Sc}\limits^{\sim }=4,$ ${C}_{g}=6.18,$ $k=-5\times {10}^{-3},$ ${n}_{1}=2,$ ${n}_{2}=3,$ ${h}_{x}=1/40,$ $\tau =1/40.$ |
Figure 4. Comparison of $u$ and $C$ with different time step size $\tau $ with $t=1,$ $\alpha =0.8,$ $\gamma =0.7,$ $\mathop{R{\rm{e}}}\limits^{\sim }=4,$ $\mathop{S{\rm{c}}}\limits^{\sim }=4,$ ${C}_{g}=6.18,$ $k=-5\times {10}^{-3},$ ${n}_{1}=2,$ ${n}_{2}=3,$ ${h}_{x}=1/40,$ ${h}_{y}=1/50.$ |
By considering the accuracy of solutions and time saving, the grid size is determined as ${h}_{x}=1/40,$ ${h}_{y}=1/50$ and $\tau =1/40$ in the following computations.
4.2. Accuracy and convergence
4.2.1. Example 1
Algorithm in section 3 is employed to solve the following coupled equations with fractional derivatives:
$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial x}+\displaystyle \frac{\partial v}{\partial y}=0,\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial u}{\partial t}=-u\displaystyle \frac{\partial u}{\partial x}-v\displaystyle \frac{\partial u}{\partial y}+{F}_{\alpha }\displaystyle \frac{{\partial }^{\alpha +1}u}{\partial {y}^{\alpha +1}}+f(x,y,t),\\ 0\lt \alpha \lt 1,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial C}{\partial t}=-u\displaystyle \frac{\partial C}{\partial x}-v\displaystyle \frac{\partial C}{\partial y}+{L}_{\gamma }\displaystyle \frac{{\partial }^{\gamma +1}C}{\partial {y}^{\gamma +1}}+g(x,y,t),\\ 0\lt \gamma \lt 1,\end{array}\end{eqnarray}$
with boundary conditions: $\begin{eqnarray}\begin{array}{l}u(x,y,0)=\displaystyle \frac{4}{3}K{x}^{\displaystyle \frac{3}{2}}y,\,v(x,y,0)=K{x}^{\displaystyle \frac{1}{2}}(1+y)\\ \,\times (1-y),C(x,y,0)=-\displaystyle \frac{\tilde{B}}{{\rm{\Gamma }}(\gamma +1)}K{x}^{\displaystyle \frac{1}{2}}{y}^{\gamma },\\ 0\leqslant x\leqslant 1,0\leqslant y\leqslant 1,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u(0,y,t)=v(0,y,t)=C(0,y,t)=0,0\leqslant y\leqslant 1,\\ 0\lt t\leqslant T,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u(x,0,t)=0,\,v(x,0,t)=K{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{1}{2}},\,\displaystyle \frac{{\partial }^{\gamma }}{\partial {y}^{\gamma }}C(x,0,t)\\ \,=\,-\tilde{B}K{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{1}{2}},\,0\leqslant x\leqslant 1,\,0\leqslant t\leqslant T,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}u(x,1,t)=\displaystyle \frac{4}{3}K{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{3}{2}},\,v(x,1,t)=0,\,C(x,1,t)\\ \,=\,-\displaystyle \frac{\tilde{B}}{{\rm{\Gamma }}(\gamma +1)}K{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{1}{2}},\,0\leqslant x\leqslant 1,\,0\leqslant t\leqslant T,\end{array}\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{l}f(x,y,t)=-\displaystyle \frac{4}{3}sK{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{3}{2}}y+\displaystyle \frac{8}{3}{K}^{2}{{\rm{e}}}^{-2st}{x}^{2}{y}^{2}\\ \,+\,\displaystyle \frac{4}{3}{K}^{2}{{\rm{e}}}^{-2st}{x}^{2}(1-{y}^{2})-{F}_{\alpha }\displaystyle \frac{4}{3}K{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{3}{2}}q(y),\end{array}\end{eqnarray}$
$\begin{eqnarray}q(y)=\displaystyle \frac{{\rm{\Gamma }}(2){y}^{-\alpha }}{{\rm{\Gamma }}(1-\alpha )},\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}g(x,y,t)=\displaystyle \frac{\tilde{B}K{{\rm{e}}}^{-st}}{{\rm{\Gamma }}(\gamma +1)}\left(s{x}^{\displaystyle \frac{1}{2}}{y}^{\gamma }-\displaystyle \frac{2K}{3}{{\rm{e}}}^{-st}x{y}^{\gamma +1}\right.\\ \,-\,\left.\gamma K{{\rm{e}}}^{-st}x({y}^{\gamma -1}-{y}^{\gamma +1})\Space{0ex}{2.8ex}{0ex}\right)+{L}_{\gamma }\tilde{B}K{{\rm{e}}}^{-st}{x}^{\displaystyle \frac{1}{2}}{y}^{-1}.\end{array}\end{eqnarray}$
The exact solution is $u(x,y,t)=\tfrac{4}{3}K{{\rm{e}}}^{-st}{x}^{\tfrac{3}{2}}y,$ $v(x,y,t)\,=K{{\rm{e}}}^{-st}{x}^{\tfrac{1}{2}}(1-y)(1+y)$ and $C(x,y,t)\,=-\tfrac{\tilde{B}K}{{\rm{\Gamma }}(\gamma +1)}{{\rm{e}}}^{-st}{x}^{\tfrac{1}{2}}{y}^{\gamma }.$
In figure 5, the calculated results are compared with known accurate solutions. The results show good consistency, and it proves the accuracy of the numerical algorithm. Therefore, we consider that the numerical algorithm herein is reliable and can be extended to other fractional-order boundary layer models with non-uniform infiltration.
Figure 5. Numerical solution and exact solution of $u,$ $v$ and $C$ with ${h}_{x}=1/40,$ ${h}_{y}=1/50$ and $\tau =1/40.$ |
Let ${u}_{i,j}^{k}$ and ${C}_{i,j}^{k}$ be the numerical solutions of the finite difference algorithm, and $u({x}_{i},{y}_{j},{t}_{k})$ and $C({x}_{i},{y}_{j},{t}_{k})$ be the exact values of the problems (25 ), (26 ). Based on the discussion in section 3 , let
$\begin{eqnarray}{e}_{i,j}^{k}={u}_{i,j}^{k}-u({x}_{i},{y}_{j},{t}_{k}),\end{eqnarray}$
$\begin{eqnarray}{\varepsilon }_{i,j}^{k}={C}_{i,j}^{k}-C({x}_{i},{y}_{j},{t}_{k}),\end{eqnarray}$
then $\begin{eqnarray}\left|{e}_{i,j}^{k}\right|\leqslant A(\tau +{h}_{x}+{h}_{y}),\end{eqnarray}$
$\begin{eqnarray}\left|{\varepsilon }_{i,j}^{k}\right|\leqslant B(\tau +{h}_{x}+{h}_{y}),\end{eqnarray}$
where $i=1,2,\,\ldots ,\,{M}_{x};$ $j=1,2,\,\ldots ,\,{M}_{y};$ $k=1,2,\,\ldots ,\,N;$ $A$ and $B$ are positive constants.Table 1. The numerical errors using the difference approximation ( |
| ${h}_{x}$ | ${h}_{y}$ | $\tau $ | ${\parallel {e}^{k}\parallel }_{\max }$ |
|---|---|---|---|
| $1/40$ | $1/40$ | $1/40$ | 2.1068e-04 |
| $1/80$ | $1/80$ | $1/80$ | 1.0726e-04 |
| $1/160$ | $1/160$ | $1/160$ | 5.3961e-05 |
| $1/200$ | $1/200$ | $1/200$ | 4.3193e-05 |
Table 2. The numerical errors using the difference approximation ( |
| ${h}_{x}$ | ${h}_{y}$ | $\tau $ | ${\parallel {\varepsilon }^{k}\parallel }_{\max }$ |
|---|---|---|---|
| $1/40$ | $1/40$ | $1/40$ | 5.3991e-03 |
| $1/80$ | $1/80$ | $1/80$ | 2.9427e-03 |
| $1/160$ | $1/160$ | $1/160$ | 1.5850e-03 |
| $1/200$ | $1/200$ | $1/200$ | 1.2975e-03 |
Let41 ) and (43 ) satisfy equations (65 ), (66 ).
$\begin{eqnarray}{\parallel {e}^{k}\parallel }_{\max }={\max }_{1\leqslant i\leqslant {M}_{x}-1,1\leqslant j\leqslant {M}_{y}-1}\left|{e}_{i,j}^{k}\right|,\end{eqnarray}$
$\begin{eqnarray}{\parallel {\varepsilon }^{k}\parallel }_{\max }={\max }_{1\leqslant i\leqslant {M}_{x}-1,1\leqslant j\leqslant {M}_{y}-1}\left|{\varepsilon }_{i,j}^{k}\right|.\end{eqnarray}$
From tables 1–2, we can conclude that the errors of the difference scheme (4.3. Comparison with other method
4.3.1. Example 2
$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial t}=-\displaystyle \frac{\partial u}{\partial x}-\displaystyle \frac{\partial u}{\partial y}+{F}_{\alpha }\displaystyle \frac{{\partial }^{\alpha }u}{\partial {y}^{\alpha }}+f(x,y,t),\,0\lt \alpha \lt 1,\end{eqnarray}$
with boundary conditions: $\begin{eqnarray}u(x,y,0)=0,\,0\leqslant x\leqslant 1,\,0\leqslant y\leqslant 1,\end{eqnarray}$
$\begin{eqnarray}u(0,y,t)=u(1,y,t)=0,\,0\leqslant y\leqslant 1,\,0\leqslant t\leqslant T,\end{eqnarray}$
$\begin{eqnarray}u(x,0,t)=u(x,1,t)=0,\,0\leqslant x\leqslant 1,\,0\leqslant t\leqslant T,\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{l}f(x,y,t)=2tx(x{(1-x)}^{2}+t(1-3x+2{x}^{2})){y}^{2}\\ \,\times \,{(1-y)}^{2}+{t}^{2}{x}^{2}{(1-x)}^{2}\left(2y(1-3y+2{y}^{2})\right.\\ \,-\,\left.{F}_{\alpha }q(y)\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}q(y)=\displaystyle \frac{{\rm{\Gamma }}(3){y}^{2-\alpha }}{{\rm{\Gamma }}(3-\alpha )}-2\displaystyle \frac{{\rm{\Gamma }}(4){y}^{3-\alpha }}{{\rm{\Gamma }}(4-\alpha )}+\displaystyle \frac{{\rm{\Gamma }}(5){y}^{4-\alpha }}{{\rm{\Gamma }}(5-\alpha )}.\end{eqnarray}$
The exact solution is $u(x,y,t)={t}^{2}{x}^{2}{(1-x)}^{2}{y}^{2}{(1-y)}^{2}.$The comparisons of the numerical results and the CPU time between the Grünwald–Letnikov approximation scheme adopted in this research and L1 approximation scheme with different fractional derivative $\alpha ,$ parameter ${F}_{\alpha }$ and mesh densities are shown in table 3. The advantages and disadvantages are quantified by these comparisons. The results of these two numerical methods are relatively consistent. It is worth noting that when the mesh generation is sparse, such as ${h}_{x}={h}_{y}=\tau =1/40$ and $1/80,$ the calculation time of L1 approximation algorithm is shorter. However, with the grid encryption, such as ${h}_{x}={h}_{y}\,=\tau =1/160,$ the calculation time of Grünwald–Letnikov approximation is obviously shorter than that of L1 algorithm. Furthermore, when Grünwald–Letnikov algorithm is employed to solve the fractional derivative equation defined by Riemann–Liouville fractional derivative, the relationship between Riemann–Liouville and Grünwald–Letnikov derivative is used to discretize the fractional derivative by the shifted Grünwald–Letnikov formula. However, L1 algorithm is always employed to discretize Caputo derivative when Riemann–Liouville derivative is transformed into Caputo type. Sun et al [57] has been pointed out that, when the fractional derivative $\alpha $ is between 1 and 2, the initial condition $u(x,y,0)=\phi (x,y)$ is not sufficient, but another initial condition $\tfrac{\partial u(x,y,0)}{\partial t}=\varphi (x,y)$ is necessary for later computation. Above all, the Grünwald–Letnikov approximation scheme is more efficient and convenient to deal with the spatial fractional derivative term with Riemann–Liouville definition.
Table 3. Comparisons between G-L algorithm and L1 algorithm with $t=1.$ |
| $(\alpha ,{F}_{\alpha })$ | ${h}_{x}={h}_{y}=\tau $ | $u$ of G-L algorithm | $u$ of L1 algorithm | CPU time of G-L algorithm | CPU time of L1 algorithm |
|---|---|---|---|---|---|
| $(0.9,0.2)$ | $1/40$ | 3.7089e-03 | 3.7380e-03 | 4.143 650 | 4.017 354 |
| $1/80$ | 3.8034e-03 | 3.8185e-03 | 52.166 716 | 51.791 007 | |
| $1/160$ | 3.8535e-03 | 3.8612e-03 | 923.833 059 | 975.019 377 | |
| $(0.8,0.4)$ | $1/40$ | 3.6876e-03 | 3.7343e-03 | 4.481 455 | 3.984 641 |
| $1/80$ | 3.7917e-03 | 3.8156e-03 | 54.090 153 | 51.002 836 | |
| $1/160$ | 3.8474e-03 | 3.8593e-03 | 928.049 600 | 954.880 134 | |
| $(0.7,0.6)$ | $1/40$ | 3.6731e-03 | 3.7271e-03 | 3.992 131 | 3.975 566 |
| $1/80$ | 3.7838e-03 | 3.8110e-03 | 52.689 535 | 50.739 646 | |
| $1/160$ | 3.8431e-03 | 3.8566e-03 | 983.034 463 | 1006.366 636 |
5. Results and discussion
In this section, the momentum and mass transfer process of viscoelastic liquid food in a permeable plate are discussed, and the influences of the time, fractional-order constitutive models, generalized Reynolds number, generalized Schmidt number and the permeability parameter on the velocity and concentration fields are considered. The physical meaning of the rate of penetration ${v}_{w}(x,t)=k{x}^{{n}_{1}}{t}^{{n}_{2}}$ can be explained according to the values of ${n}_{1}$ and ${n}_{2}.$ For instance, when ${n}_{1}=2,$ ${n}_{2}=3,$ it is assumed that the permeability pores at the entrance of the channel are small, and then gradually enlarge with the increase of the distance along the $x$ direction. Thus, the permeability enhances with the increase of $x.$ On the other hand, with the passage of time, the amount of fluid flowing into the filter channel gradually enlarges, resulting in a bigger penetration rate of the fluid.
The impact of time on the velocity and concentration distribution in the boundary layer is studied in figure 6, revealing the three-dimensional distribution of velocity and concentration. Figures 6(a) and (b) are the three-dimensional distribution of velocity $u$ with respect to time $t$ in $x$ and $y$ direction, respectively. Figures 6(c) and (d) are the three-dimensional distribution of concentration $C$ with respect to time $t$ in $x$ and $y$ direction, respectively. Initially, the change in the fluid velocity and concentration is more intense. With the passage of time, the fluid velocity and concentration gradually decrease and finally tend to be flat. The fluids flow transitions from an unsteady state to a steady one and finally reaches a stable state.
Figure 6. Comparison of $u$ and $C$ with different $t$ for $\alpha =0.8,$ $\gamma =0.7,$ $\mathop{R{\rm{e}}}\limits^{\sim }=4,$ $\mathop{Sc}\limits^{\sim }=4,$ ${C}_{g}=6.18,$ $k=-5\times {10}^{-5},$ ${n}_{1}=2,$ ${n}_{2}=3:$ (a) $y=0.2;$ (b) $x=0.4;$ (c) $y=0.2;$ (d) $x=0.4.$ |
Figures 7 and 8 demonstrate the effect of fractional-order constitutive models, i.e. models (4 ) and (6 ), the fractional derivatives ($0.8\leqslant \alpha \leqslant 1.0,$ $0.7\leqslant \gamma \leqslant 1.0$) on the velocity and concentration boundary layers. In figure 7, the increase in fractional-order $\alpha $ results in a decrease in horizontal velocity $u.$ This is consistent with [27, 28]. The higher the fractional order $\alpha $ is; the thicker the momentum boundary layer is; the greater the flow resistance of the fluid is, which in turn decreases the fluid velocity. Figure 8 describes the relation between concentration and fractional derivative $\gamma .$ Overall, the increase in fractional-order $\gamma $ results in a decrease in concentration $C,$ because anomalous diffusion weakens the mass transfer. In short, the increase in fractional derivative has a significant inhibitory effect on mass transfer. Furthermore, figure 8(a) compares the concentration distribution of the fluid with varying fractional derivatives $\gamma $ at $y=0.1$ and $y=0.5,$ respectively. The concentration is lower at points closer to the edge of the plate than that at the middle of the boundary layer; this is due to the influence of wall permeability.
Figure 7. Comparison of $u$ with varying $\alpha $ for $t=1,$ $\tilde{R{\rm{e}}}=4,$ $k=-5\times {10}^{-4},$ ${n}_{1}=2,$ ${n}_{2}=3.$ |
Figure 8. Comparison of $C$ with different $\gamma $ for $t=1,$ $\alpha =0.8,$ $\mathop{R{\rm{e}}}\limits^{\sim }=4,$ $\mathop{Sc}\limits^{\sim }=4,$ ${C}_{g}=6.18,$ $k=-5\times {10}^{-{\rm{4}}},$ ${n}_{1}=2,$ ${n}_{2}=3.$ |
Figure 9 is a distribution diagram used to examine the influence of generalized Reynolds number on velocity fields. It is noted that the generalized Reynolds number in the fractional differential equation depicting laminar flow is not as large as that in the differential equation of integer-order. In [26], Sun et al mentioned that fractional generalized Reynolds number is different from integer-order Reynolds number. Pan et al [27] examined the influence of generalized Reynolds number on velocity profile and the range of Reynolds number is limited to $1\leqslant \mathop{R{\rm{e}}}\limits^{\sim }\leqslant 2.$ Tabi et al [58] showed that when $\mathop{R{\rm{e}}}\limits^{\sim }=3,$ $4,$ $5,$ the flow velocity of fluid increased significantly with the increase of $\mathop{R{\rm{e}}}\limits^{\sim }.$ In addition, the range of Reynolds numbers were chosen as $\mathop{R{\rm{e}}}\limits^{\sim }=2,$ $3,$ $6,$ $10$ in the investigation on blood flow characteristics, considered by Maiti et al [59]. Adopting similar strategy as above literatures, the range of generalized Reynolds number in this paper is $1\leqslant \mathop{R{\rm{e}}}\limits^{\sim }\leqslant 10.$ In figure 9(a), as $\mathop{R{\rm{e}}}\limits^{\sim }$ increases, the flow viscous coefficient ${\tilde{\mu }}_{\alpha }$ decreases, resulting in the reduction of thickness in the velocity boundary layer and an increase in fluid velocity. In figure 9(b), the velocity distributions of the fluid at $y=0.1$ and $y=0.5$ with varying generalized Reynolds number are compared. Due to the influence of wall permeability, the velocity near the plate is much lower than the velocity in the middle of the boundary layer.
Figure 9. Comparison of $u$ with varying $\mathop{R{\rm{e}}}\limits^{\sim }$ for $t=1,$ $\alpha =0.8,$ $k=-5\times {10}^{-4},$ ${n}_{1}=2,$ ${n}_{2}=3.$ |
Figure 10 observes the impact of generalized Schmidt number on the concentration distribution. The ratio of the kinematic viscosity to the molecular diffusion is defined as the generalized Schmidt number. When $\mathop{Sc}\limits^{\sim }$ increases, the molecular diffusion coefficient ${\tilde{D}}_{\gamma }$ reduces, resulting in a decrease in concentration boundary layer thickness and an increase in fluid concentration, which is confirmed by figure 10. Figure 10(a) shows an interesting phenomenon that there is a small reversal of fluid concentration near the wall. The concentration first decreases a little and then gradually increases along the $y$ direction. When the plate is in the state of seepage, part of the fluid will permeate through the permeable plate, resulting in the decrease in fluid concentration. However, with an increase in $y,$ the fluid flow is affected strongly by the main flow, bringing about a gradual increase in fluid concentration, which is close to the dimensionless concentration 1 in the end. In figure 10(b), the concentrations along the $x$ direction with varying Schmidt numbers are compared at $y=0.1$ and $y=0.5.$ The fluid concentration drops sharply on the inlet, then rises slightly and finally tends to be gentle with an increase in $x.$ At the inlet, the fluid is suddenly constrained by the plate surface; thus, the concentration changes from $1$ to the boundary layer concentration. With the enhancement in $x,$ the permeability becomes stronger, which causes the fluid in the boundary layer appreciably affected by the mainstream. Consequently, the concentration of the boundary layer increases slightly.
Figure 10. Comparison of $C$ with different $\mathop{Sc}\limits^{\sim }$ for $t=1,$ $\alpha =0.8,$ $\gamma =0.7,$ $\mathop{R{\rm{e}}}\limits^{\sim }=4,$ ${C}_{g}=6.18,$ $k=-5\times {10}^{-{\rm{4}}},$ ${n}_{1}=2,$ ${n}_{2}=3.$ |
To examine the impact of permeability on the mass transfer, figure 11 shows the concentration profiles along the $x$ direction by changing the value of the permeability parameter $k.$ In equation (24 ), ${v}_{w}(x,t)$ is written as a specific function. $k\lt 0$ corresponds to the exudation, and $k\gt 0$ corresponds to the injection. With an increase in $x,$ the fluid concentration begins to sharply decrease and then gradually be gentle. This is because near the entrance of the channel, the fluid in the boundary layer is suddenly restrained by the plate; thus, the concentration drops sharply from the dimensionless concentration 1 of the main stream to a concentration of approximately 0.3. Subsequently, the concentration decreased/increased slowly under the influence of the plate permeability. With the increase in $x,$ the permeability velocity ${v}_{w}(x,t)$ increases gradually, leading to an increase in the quantity of fluid exuded/injected into the permeable plate. When $x$ tends to be the dimensionless number 1, the influence of permeability characteristics on the concentration of the boundary layer reaches the maximum. In addition, when the absolute value of $k$ is very small, for example, $k=5\times {10}^{-7},$ the influence of the positive and negative permeability parameters on mass transfer is inconspicuous.
Figure 11. Comparison of $C$ with different $k$ for $y=0.06,$ $t=3,$ $\alpha =0.8,$ $\gamma =0.7,$ $\mathop{R{\rm{e}}}\limits^{\sim }=4,$ ${n}_{1}=2,$ ${n}_{2}=3.$ |
As mentioned earlier, anomalous diffusion is very common in food and biomaterial processing. The implicit finite difference method adopted in this research is employed to solve a similar fractional diffusion equation to stimulate the diffusion of nisin in 3% agarose gel column in 6.11 days with a constant experimental temperature of 9.9 °C. The diffusion of nisin in agarose gels was also characterized by the modified Fick's law ${J}_{y}={\tilde{D}}_{\gamma }\tfrac{{\partial }^{\gamma }C}{\partial {y}^{\gamma }}.$ A comparison between the experimental data from [60] and the numerical results obtained by us has been revealed in figure 12. According to [60], the corresponding diffusion coefficient is $3.52\times {10}^{-11}$ $({{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}),$ and the parameters used in the numerical calculation are the same as those in the literature. Figure 12 shows our numerical results are in a good agreement with the experimental data, which further proves the reliability of the model and method in this paper.
Figure 12. Comparison of experimental data from [60] with our results to test the method and the model. |
6. Conclusions
In this study, the boundary layer flow of liquid food on a non-uniform permeable plate is investigated. A fractional constitutive model is introduced to describe the viscoelastic fluid and to generalize Fick's law. The fluid on the plate permeates at an uneven rate, as the pores may be blocked by the settled impurities in the fluid, or the liquid food may spontaneously form a deposition of a gel layer. Considering the steady-state mass balance during ultrafiltration on the plate surface, a fractional-order concentration boundary condition is established. The finite difference algorithm is used to solve the problem, and its accuracy and grid independence are verified, proving that the numerical method can efficiently solve similar spatial fractional differential equations. With the aid of the numerical algorithm, some interesting theoretical analysis for the liquid filtration, which is essential in food processing, is observed. The effects of time, fractional constitutive models, generalized Reynolds number, generalized Schmidt number, and permeability parameter on the velocity and concentration boundary layers are revealed. The main conclusions are as follows:
| 1. An increase in fractional-order $\alpha $ brings about a decrease in the horizontal velocity. | |
| 2. Anomalous diffusion weakens the mass transfer; therefore, the concentration decreases by increasing the fractional derivative $\gamma .$ | |
| 3. With an increase in the generalized Reynolds number, the boundary layer thickness reduces, while the fluid velocity increases. | |
| 4. When the generalized Schmidt number rises, the thickness of the concentration boundary layer drops, and the fluid concentration increases. | |
| 5. It is also worth noting that the concentration suddenly decreases slightly near the wall and then gradually increases with an increase in $y.$ |