On the process of filtration of fractional viscoelastic liquid food

Yahui Meng; Botong Li; Xinhui Si; Xuehui Chen; Fawang Liu

doi:10.1088/1572-9494/abdf45

PDF(1474 KB)

Communications in Theoretical Physics ›› 2021, Vol. 73 ›› Issue (4) : 45004. DOI: 10.1088/1572-9494/abdf45

Mathematical Physics

On the process of filtration of fractional viscoelastic liquid food

Author information +

History +

Abstract

In the process of filtration, fluid impurities precipitate/accumulate; this results in an uneven inner wall of the filter, consequently leading to non-uniform suction/injection. The Riemannian–Liouville fractional derivative model is used to investigate viscoelastic incompressible liquid food flowing through a permeable plate and to generalize Fick's law. Moreover, we consider steady-state mass balance during ultrafiltration on a plate surface, and a fractional-order concentration boundary condition is established, thereby rendering the problem real and complex. The governing equation is numerically solved using the finite difference algorithm. The effects of the fractional constitutive models, generalized Reynolds number, generalized Schmidt number, and permeability parameter on the velocity and concentration fields are compared. The results show that an increase in fractional-order α in the momentum equation leads to a decrease in the horizontal velocity. Anomalous diffusion described by the fractional derivative model weakens the mass transfer; therefore, the concentration decreases with increasing fractional derivative γ in the concentration equation.

Key words

non-uniform permeable surface / viscoelastic fluid / fractional calculus / implicit numerical method

Cite this article

EndNote

Ris (Procite)

Bibtex

Download Citations

Yahui Meng, Botong Li, Xinhui Si, et al. On the process of filtration of fractional viscoelastic liquid food[J]. Communications in Theoretical Physics, 2021, 73(4): 45004 https://doi.org/10.1088/1572-9494/abdf45

Nomenclature

$C$	fluid concentration
$C_{\infty}$	concentration at infinity
$C_{g}$	gel layer concentration
${\tilde{D}}_{γ}$	generalized diffusion coefficient
$J_{y}$	mass transfer flux
$k$	permeability parameter
$L$	characteristic length
$n_{1}, n_{2}$	constants
$\tilde{R e}$	generalized Reynolds number
$\tilde{S c}$	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
$α$	fractional spatial derivative of velocity
$γ$	fractional spatial derivative of concentration
${\tilde{μ}}_{α}$	generalized dynamic viscosity
$ρ$	fluid density
$τ_{y x}$	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 Fe₃O₄-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} + {D \frac{\partial C}{\partial y} |}_{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} + \tilde{D_{γ}} \frac{\partial^{γ} C}{\partial y^{γ}} (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.

Full size|PPT slide

The governing equations are given as follows.

Continuity equation:

\begin{array}{rcl} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0. \end{array}

(1)

Momentum and concentration equations:

\begin{array}{rcl} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \frac{1}{ρ} \frac{\partial τ_{y x}}{\partial y}, \end{array}

(2)

\begin{array}{rcl} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = \frac{\partial J_{y}}{\partial y} . \end{array}

(3)

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{array}{rcl} τ_{y x} = {\tilde{μ}}_{α} \frac{\partial^{α} u}{\partial y^{α}}, \end{array}

(4)

where

{\tilde{μ}}_{α}

denotes the generalized dynamic viscosity, and

α

(0 < α < 1)

is the space-fractional derivative. The Riemann–Liouville fractional derivative is used [49]:

\begin{array}{rcl} \frac{\partial^{α} u}{\partial y^{α}} = \frac{1}{Γ (1 - α)} \frac{\partial}{\partial y} \int_{0}^{y} {(y - ξ)}^{- α} u (x, ξ, t) d ξ, \end{array}

(5)

where

Γ (\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⁻¹ 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{array}{rcl} J_{y} = {\tilde{D}}_{γ} \frac{\partial^{γ} C}{\partial y^{γ}}, \end{array}

(6)

where

{\tilde{D}}_{γ}

denotes the generalized diffusion coefficient, and

γ

(0 < γ < 1)

is the fractional spatial derivative. The Riemann–Liouville fractional derivative is also adopted as in [49]:

\begin{array}{rcl} \frac{\partial^{γ} C}{\partial y^{γ}} = \frac{1}{Γ (1 - γ)} \frac{\partial}{\partial y} \int_{0}^{y} {(y - ξ)}^{- γ} C (x, ξ, t) d ξ . \end{array}

(7)

By using the fractional-order relations of (4) and (6), the governing equations (2) and (3) become

\begin{array}{rcl} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \frac{{\tilde{μ}}_{α}}{ρ} \frac{\partial}{\partial y} (\frac{\partial^{α} u}{\partial y^{α}}), \end{array}

(8)

\begin{array}{rcl} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = {\tilde{D}}_{γ} \frac{\partial}{\partial y} (\frac{\partial^{γ} C}{\partial y^{γ}}) . \end{array}

(9)

The initial and boundary conditions are as follows:

\begin{array}{rcl} \begin{array}{l} u (x, y, 0) = u_{\infty}, v (x, y, 0) = 0, C (x, y, 0) = C_{\infty}, \\ x ⩾ 0, y > 0, \end{array} \end{array}

(10)

\begin{array}{rcl} u (x, 0, t) = 0, v (x, 0, t) = v_{w} (x, t), x > 0, t > 0, \end{array}

(11)

\begin{array}{rcl} v_{w} (x, t) C_{g} + \tilde{D_{γ}} \frac{\partial^{γ} C}{\partial y^{γ}} (x, 0, t) = 0, x > 0, t > 0, \end{array}

(12)

\begin{array}{rcl} \begin{array}{l} u (x, + \infty, t) = u_{\infty}, v (x, + \infty, t) = 0, \\ C (x, + \infty, t) = C_{\infty}, x ⩾ 0, t > 0, \end{array} \end{array}

(13)

\begin{array}{rcl} \begin{array}{l} u (0, y, t) = u_{\infty}, v (0, y, t) = 0, C (0, y, t) = C_{\infty}, \\ y > 0, t ⩾ 0, \end{array} \end{array}

(14)

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 < 0

corresponds to suction, and

k > 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} + {D \frac{\partial C}{\partial y} |}_{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 \frac{\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}}_{γ} \frac{\partial^{γ} C}{\partial y^{γ}}

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 (12).

C_{g}

is the concentration of the gel layer.

The following dimensionless variables are adopted:

\begin{array}{rcl} \begin{array}{l} x * = \frac{x}{L}, y * = \frac{y}{L}, u * = \frac{u}{u_{\infty}}, v * = \frac{v}{u_{\infty}}, v_{w}^{*} = \frac{v_{w}}{u_{\infty}}, \\ C * = \frac{C}{C_{\infty}}, C_{g}^{*} = \frac{C_{g}}{C_{\infty}}, t * = \frac{t u_{\infty}}{L}, k * = \frac{L^{n_{1} + n_{2}}}{u_{\infty}^{n_{2} + 1}} k, \end{array} \end{array}

(15)

where characteristic length is denoted by

L .

The dimensionless governing equations become

\begin{array}{rcl} \frac{\partial u *}{\partial x *} + \frac{\partial v *}{\partial y *} = 0, \end{array}

(16)

\begin{array}{rcl} \frac{\partial u *}{\partial t *} + u * \frac{\partial u *}{\partial x *} + v * \frac{\partial u *}{\partial y *} = \frac{1}{\tilde{R e}} \frac{\partial}{\partial y *} (\frac{\partial^{α} u *}{\partial y^{* α}}), \end{array}

(17)

\begin{array}{rcl} \frac{\partial C *}{\partial t *} + u * \frac{\partial C *}{\partial x *} + v * \frac{\partial C *}{\partial y *} = \frac{1}{\tilde{R e}} \frac{1}{\tilde{S c}} \frac{\partial}{\partial y *} (\frac{\partial^{γ} C *}{\partial y^{* γ}}), \end{array}

(18)

where

\tilde{R e} = \frac{ρ u_{\infty} L^{α}}{{\tilde{μ}}_{α}}

is the generalized Reynolds number, and

\tilde{S c} = \frac{{\tilde{μ}}_{α}}{ρ L^{α - γ} {\tilde{D}}_{γ}}

is the generalized Schmidt number.

The corresponding non-dimensional initial and boundary conditions are

\begin{array}{rcl} \begin{array}{l} u * (x *, y *, 0) = 1, v * (x *, y *, 0) = 0, \\ C * (x *, y *, 0) = 1, x * ⩾ 0, y * > 0, \end{array} \end{array}

(19)

\begin{array}{rcl} \begin{array}{l} u * (x *, 0, t *) = 0, v * (x *, 0, t *) = v_{w}^{*} (x *, t *), \\ x * > 0, t * > 0, \end{array} \end{array}

(20)

\begin{array}{rcl} \begin{array}{l} \tilde{R e} \tilde{S c} C_{g}^{*} v_{w}^{*} (x *, t *) + \frac{\partial^{γ} C *}{\partial y^{* γ}} (x *, 0, t *) = 0, \\ x * > 0, t * > 0, \end{array} \end{array}

(21)

\begin{array}{rcl} \begin{array}{l} u * (x *, + \infty, t *) = 1, v * (x *, + \infty, t *) = 0, \\ C * (x *, + \infty, t *) = 1, x * ⩾ 0, t * > 0, \end{array} \end{array}

(22)

\begin{array}{rcl} \begin{array}{l} u * (0, y *, t *) = 1, v * (0, y *, t *) = 0, \\ C * (0, y *, t *) = 1, y * > 0, t * ⩾ 0, \end{array} \end{array}

(23)

where

\begin{array}{rcl} v_{w}^{*} (x *, t *) = k * x^{* n_{1}} t^{* n_{2}} . \end{array}

(24)

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) as

\begin{array}{rcl} \begin{array}{l} \frac{\partial u}{\partial t} + a (x, y, t) \frac{\partial u}{\partial x} + b (x, y, t) \frac{\partial u}{\partial y} = F_{α} \frac{\partial^{α + 1} u}{\partial y^{α + 1}} \\ + f (x, y, t), \end{array} \end{array}

(25)

\begin{array}{rcl} \begin{array}{l} \frac{\partial C}{\partial t} + a (x, y, t) \frac{\partial C}{\partial x} + b (x, y, t) \frac{\partial C}{\partial y} = L_{γ} \frac{\partial^{γ + 1} C}{\partial y^{γ + 1}} \\ + g (x, y, t), \end{array} \end{array}

(26)

where

F_{α} = \frac{1}{\tilde{R e}},

L_{γ} = \frac{1}{\tilde{R e}} \frac{1}{\tilde{S c}} .

In the physical problem presented in section 2,

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{array}{rcl} \frac{\partial a (x, y, t)}{\partial x} + \frac{\partial b (x, y, t)}{\partial y} = 0. \end{array}

(27)

Define

t_{k} = k τ,

k = 0, 1, \dots, n,

x_{i} = i h_{x},

i = 0, 1, \dots, M_{x},

y_{j} = j h_{y},

j = 0, 1, \dots, M_{y},

where

τ = 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^{α + 1} u / \partial y^{α + 1}

and

\partial^{γ + 1} C / \partial y^{γ + 1}

are as follows [56]:

\begin{array}{rcl} \begin{array}{l} \frac{\partial^{α + 1}}{\partial y^{α + 1}} u (x_{i}, y_{j}, t_{k + 1}) = \frac{1}{h_{y}^{α + 1}} \sum_{m = 0}^{j + 1} ω_{m}^{α + 1} u (x_{i}, y_{j + 1 - m}, \\ t_{k + 1}) + O (h_{y}), \end{array} \end{array}

(28)

\begin{array}{rcl} \begin{array}{l} \frac{\partial^{γ + 1}}{\partial y^{γ + 1}} C (x_{i}, y_{j}, t_{k + 1}) = \frac{1}{h_{y}^{γ + 1}} \sum_{q = 0}^{j + 1} η_{q}^{γ + 1} C (x_{i}, y_{j + 1 - q}, \\ t_{k + 1}) + O (h_{y}) . \end{array} \end{array}

(29)

Taking as

\begin{array}{rcl} \begin{array}{l} ω_{0}^{α + 1} = 1, ω_{m}^{α + 1} = (1 - \frac{α + 2}{m}) ω_{m - 1}^{α + 1}, m = 1, 2, \dots, \\ η_{0}^{γ + 1} = 1, η_{q}^{γ + 1} = (1 - \frac{γ + 2}{q}) η_{q - 1}^{γ + 1}, q = 1, 2, \dots . \end{array} \end{array}

Using the backward difference method,

\partial u / \partial x

and

\partial C / \partial x

can be approximated as

\begin{array}{rcl} \begin{array}{l} \frac{\partial u (x_{i}, y_{j}, t_{k + 1})}{\partial x} = \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{array}

(30)

\begin{array}{rcl} \begin{array}{l} \frac{\partial C (x_{i}, y_{j}, t_{k + 1})}{\partial x} = \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{array}

(31)

Similarly,

\partial u / \partial y

and

\partial C / \partial y

are approximated as

\begin{array}{rcl} \begin{array}{l} \frac{\partial u (x_{i}, y_{j}, t_{k + 1})}{\partial y} = \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{array}

(32)

\begin{array}{rcl} \begin{array}{l} \frac{\partial C (x_{i}, y_{j}, t_{k + 1})}{\partial y} = \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{array}

(33)

For the first-order temporal derivative, the approximate discrete scheme is obtained by using the backward difference scheme:

\begin{array}{rcl} \begin{array}{l} \frac{\partial u (x_{i}, y_{j}, t_{k + 1})}{\partial t} = \frac{u (x_{i}, y_{j}, t_{k + 1}) - u (x_{i}, y_{j}, t_{k})}{τ} \\ + O (τ), \end{array} \end{array}

(34)

\begin{array}{rcl} \begin{array}{l} \frac{\partial C (x_{i}, y_{j}, t_{k + 1})}{\partial t} = \frac{C (x_{i}, y_{j}, t_{k + 1}) - C (x_{i}, y_{j}, t_{k})}{τ} \\ + O (τ) . \end{array} \end{array}

(35)

Hence, equations (25) and (26) become

\begin{array}{rcl} \begin{array}{l} \frac{u (x_{i}, y_{j}, t_{k + 1}) - u (x_{i}, y_{j}, t_{k})}{τ} = - a (x_{i}, y_{j}, t_{k}) \\ \times \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 \frac{u (x_{i}, y_{j}, t_{k + 1}) - u (x_{i}, y_{j - 1}, t_{k + 1})}{h_{y}} \\ + \frac{F_{α}}{h_{y}^{α + 1}} \sum_{m = 0}^{j + 1} ω_{m}^{α + 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{array}

(36)

and

\begin{array}{rcl} \begin{array}{l} \frac{C (x_{i}, y_{j}, t_{k + 1}) - C (x_{i}, y_{j}, t_{k})}{τ} = - a (x_{i}, y_{j}, t_{k + 1}) \\ \times \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}) \frac{C (x_{i}, y_{j}, t_{k + 1}) - C (x_{i}, y_{j - 1}, t_{k + 1})}{h_{y}} \\ + \frac{L_{γ}}{h_{y}^{γ + 1}} \sum_{q = 0}^{j + 1} η_{q}^{γ + 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{array}

(37)

where

\begin{array}{rcl} | R_{i, j, k + 1} | ⩽ A (τ + h_{x} + h_{y}), \end{array}

(38)

\begin{array}{rcl} | Z_{i, j, k + 1} | ⩽ B (τ + h_{x} + h_{y}) . \end{array}

(39)

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 (25) utilizing the implicit difference scheme is as follows:

\begin{array}{rcl} \begin{array}{l} \frac{u_{i, j}^{k + 1} - u_{i, j}^{k}}{τ} = - a_{i, j}^{k} \frac{u_{i, j}^{k + 1} - u_{i - 1, j}^{k + 1}}{h_{x}} - b_{i, j}^{k} \frac{u_{i, j}^{k + 1} - u_{i, j - 1}^{k + 1}}{h_{y}} \\ + \frac{F_{α}}{h_{y}^{α + 1}} \sum_{m = 0}^{j + 1} ω_{m}^{α + 1} u_{i, j + 1 - m}^{k + 1} + f_{i, j}^{k + 1}, \end{array} \end{array}

(40)

\begin{array}{rcl} \begin{array}{l} μ_{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} \sum_{m = 0, m \neq 1}^{j + 1} ω_{m}^{α + 1} \\ \times u_{i, j + 1 - m}^{k + 1} = u_{i, j}^{k} + τ f_{i, j}^{k + 1}, \end{array} \end{array}

(41)

where

r_{1} = \frac{τ}{h_{x}},

r_{2} = \frac{τ}{h_{y}},

r_{3} = \frac{τ F_{α}}{h_{y}^{α + 1}},

μ_{i, j}^{k} = 1 + r_{1} a_{i, j}^{k} + r_{2} b_{i, j}^{k} - r_{3} ω_{1}^{α + 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{array}{rcl} \begin{array}{l} \frac{C_{i, j}^{k + 1} - C_{i, j}^{k}}{τ} = - a_{i, j}^{k + 1} \frac{C_{i, j}^{k + 1} - C_{i - 1, j}^{k + 1}}{h_{x}} - b_{i, j}^{k + 1} \\ \times \frac{C_{i, j}^{k + 1} - C_{i, j - 1}^{k + 1}}{h_{y}} + \frac{L_{γ}}{h_{y}^{γ + 1}} \sum_{q = 0}^{j + 1} η_{q}^{γ + 1} C_{i, j + 1 - q}^{k + 1} + g_{i, j}^{k + 1}, \end{array} \end{array}

(42)

\begin{array}{rcl} \begin{array}{l} λ_{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} \sum_{q = 0, q \neq 1}^{j + 1} η_{q}^{γ + 1} \\ \times C_{i, j + 1 - q}^{k + 1} = C_{i, j}^{k} + τ g_{i, j}^{k + 1}, \end{array} \end{array}

(43)

where

r_{1} = \frac{τ}{h_{x}},

r_{2} = \frac{τ}{h_{y}},

r_{4} = \frac{τ L_{γ}}{h_{y}^{γ + 1}},

λ_{i, j}^{k + 1} = 1 + r_{1} a_{i, j}^{k + 1} + r_{2} b_{i, j}^{k + 1} - r_{4} η_{1}^{γ + 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^{γ} C / \partial y^{γ}

are also adopted, i.e.

\begin{array}{rcl} \begin{array}{l} \frac{\partial^{γ}}{\partial y^{γ}} C (x_{i}, y_{j}, t_{k + 1}) = \frac{1}{h_{y}^{γ}} \sum_{p = 0}^{j + 1} χ_{p}^{γ} C (x_{i}, y_{j + 1 - p}, t_{k + 1}) \\ + O (h_{y}) . \end{array} \end{array}

(44)

Taking as

χ_{0}^{γ} = 1,

χ_{p}^{γ} = (1 - \frac{γ + 1}{p}) χ_{p - 1}^{γ},

p = 1, 2, \dots .

Boundary condition (21) becomes

\begin{array}{rcl} \begin{array}{l} \frac{1}{h_{y}^{γ}} \sum_{p = 0}^{j + 1} χ_{p}^{γ} C (x_{i}, y_{j + 1 - p}, t_{k + 1}) = - \tilde{R e} \tilde{S c} C_{g} v_{w} (x_{i}, t_{k + 1}) \\ + M_{i, j, k + 1}, j = 1, \end{array} \end{array}

(45)

where

\begin{array}{rcl} | M_{i, j, k + 1} | ⩽ N (τ + h_{x} + h_{y}) . \end{array}

(46)

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:

\begin{array}{rcl} d_{1} \sum_{p = 0}^{j + 1} χ_{p}^{γ} C_{i, j + 1 - p}^{k + 1} = - \tilde{B} v_{w i}^{k + 1}, j = 1, \end{array}

(47)

\begin{array}{rcl} C_{i, 0}^{k + 1} = (- \frac{\tilde{B}}{d_{1}} v_{w i}^{k + 1} - χ_{0}^{γ} C_{i, 1}^{k + 1}) / χ_{1}^{γ}, \end{array}

(48)

where

\tilde{B} = \tilde{R e} \tilde{S c} C_{g},

d_{1} = \frac{1}{h_{y}^{γ}} .

Combining equations (43) and (48), the following equation is obtained:

\begin{array}{rcl} \begin{array}{l} λ_{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 \sum_{q = 0, q \neq 1}^{j} η_{q}^{γ + 1} C_{i, j + 1 - q}^{k + 1} - r_{4} η_{j + 1}^{γ + 1} C_{i, 0}^{k + 1} = C_{i, j}^{k} + τ g_{i, j}^{k + 1}, \end{array} \end{array}

(49)

i.e.

\begin{array}{rcl} \begin{array}{l} λ_{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} \sum_{q = 0, q \neq 1}^{j} η_{q}^{γ + 1} \\ \times C_{i, j + 1 - q}^{k + 1} + r_{4} η_{j + 1}^{γ + 1} \frac{χ_{0}^{γ}}{χ_{1}^{γ}} C_{i, 1}^{k + 1} = C_{i, j}^{k} + τ g_{i, j}^{k + 1} - θ_{i, j}^{k + 1}, \end{array} \end{array}

(50)

where

θ_{i, j}^{k + 1} = r_{4} η_{j + 1}^{γ + 1} \frac{\tilde{B}}{d_{1} χ_{1}^{γ}} v_{w i}^{k + 1} .

Note that

v (x_{i}, y_{j}, t_{k + 1})

is obtained from equation (27):

\begin{array}{rcl} \begin{array}{l} \frac{u (x_{i}, y_{j}, t_{k + 1}) - u (x_{i - 1}, y_{j}, t_{k + 1})}{h_{x}} + O (h_{x}) \\ = - \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{array}

(51)

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{array}{rcl} v_{i, j}^{k + 1} = \frac{h_{y}}{h_{x}} (u_{i - 1, j}^{k + 1} - u_{i, j}^{k + 1}) + v_{i, j - 1}^{k + 1} . \end{array}

(52)

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

τ

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,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{y} = 1 / 50,$ $τ = 1 / 50.$

Full size|PPT slide

Figure 3. Comparison of $u$ and $C$ with different space step size $h_{y}$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $τ = 1 / 40.$

Full size|PPT slide

Figure 4. Comparison of $u$ and $C$ with different time step size $τ$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $h_{y} = 1 / 50.$

Full size|PPT slide

By considering the accuracy of solutions and time saving, the grid size is determined as

h_{x} = 1 / 40,

h_{y} = 1 / 50

and

τ = 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{array}{rcl} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0, \end{array}

(53)

\begin{array}{rcl} \begin{array}{l} \frac{\partial u}{\partial t} = - u \frac{\partial u}{\partial x} - v \frac{\partial u}{\partial y} + F_{α} \frac{\partial^{α + 1} u}{\partial y^{α + 1}} + f (x, y, t), \\ 0 < α < 1, \end{array} \end{array}

(54)

\begin{array}{rcl} \begin{array}{l} \frac{\partial C}{\partial t} = - u \frac{\partial C}{\partial x} - v \frac{\partial C}{\partial y} + L_{γ} \frac{\partial^{γ + 1} C}{\partial y^{γ + 1}} + g (x, y, t), \\ 0 < γ < 1, \end{array} \end{array}

(55)

with boundary conditions:

\begin{array}{rcl} \begin{array}{l} u (x, y, 0) = \frac{4}{3} K x^{\frac{3}{2}} y, v (x, y, 0) = K x^{\frac{1}{2}} (1 + y) \\ \times (1 - y), C (x, y, 0) = - \frac{\tilde{B}}{Γ (γ + 1)} K x^{\frac{1}{2}} y^{γ}, \\ 0 ⩽ x ⩽ 1, 0 ⩽ y ⩽ 1, \end{array} \end{array}

(56)

\begin{array}{rcl} \begin{array}{l} u (0, y, t) = v (0, y, t) = C (0, y, t) = 0, 0 ⩽ y ⩽ 1, \\ 0 < t ⩽ T, \end{array} \end{array}

(57)

\begin{array}{rcl} \begin{array}{l} u (x, 0, t) = 0, v (x, 0, t) = K e^{- s t} x^{\frac{1}{2}}, \frac{\partial^{γ}}{\partial y^{γ}} C (x, 0, t) \\ = - \tilde{B} K e^{- s t} x^{\frac{1}{2}}, 0 ⩽ x ⩽ 1, 0 ⩽ t ⩽ T, \end{array} \end{array}

(58)

\begin{array}{rcl} \begin{array}{l} u (x, 1, t) = \frac{4}{3} K e^{- s t} x^{\frac{3}{2}}, v (x, 1, t) = 0, C (x, 1, t) \\ = - \frac{\tilde{B}}{Γ (γ + 1)} K e^{- s t} x^{\frac{1}{2}}, 0 ⩽ x ⩽ 1, 0 ⩽ t ⩽ T, \end{array} \end{array}

(59)

where

\begin{array}{rcl} \begin{array}{l} f (x, y, t) = - \frac{4}{3} s K e^{- s t} x^{\frac{3}{2}} y + \frac{8}{3} K^{2} e^{- 2 s t} x^{2} y^{2} \\ + \frac{4}{3} K^{2} e^{- 2 s t} x^{2} (1 - y^{2}) - F_{α} \frac{4}{3} K e^{- s t} x^{\frac{3}{2}} q (y), \end{array} \end{array}

(60)

\begin{array}{rcl} q (y) = \frac{Γ (2) y^{- α}}{Γ (1 - α)}, \end{array}

(61)

\begin{array}{rcl} \begin{array}{l} g (x, y, t) = \frac{\tilde{B} K e^{- s t}}{Γ (γ + 1)} (s x^{\frac{1}{2}} y^{γ} - \frac{2 K}{3} e^{- s t} x y^{γ + 1} \\ - γ K e^{- s t} x (y^{γ - 1} - y^{γ + 1})) + L_{γ} \tilde{B} K e^{- s t} x^{\frac{1}{2}} y^{- 1} . \end{array} \end{array}

(62)

The exact solution is

u (x, y, t) = \frac{4}{3} K e^{- s t} x^{\frac{3}{2}} y,

v (x, y, t) = K e^{- s t} x^{\frac{1}{2}} (1 - y) (1 + y)

and

C (x, y, t) = - \frac{\tilde{B} K}{Γ (γ + 1)} e^{- s t} x^{\frac{1}{2}} y^{γ} .

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 $τ = 1 / 40.$

Full size|PPT slide

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{array}{rcl} e_{i, j}^{k} = u_{i, j}^{k} - u (x_{i}, y_{j}, t_{k}), \end{array}

(63)

\begin{array}{rcl} ε_{i, j}^{k} = C_{i, j}^{k} - C (x_{i}, y_{j}, t_{k}), \end{array}

(64)

then

\begin{array}{rcl} | e_{i, j}^{k} | ⩽ A (τ + h_{x} + h_{y}), \end{array}

(65)

\begin{array}{rcl} | ε_{i, j}^{k} | ⩽ B (τ + h_{x} + h_{y}), \end{array}

(66)

where

i = 1, 2, \dots, M_{x};

j = 1, 2, \dots, M_{y};

k = 1, 2, \dots, N;

A

and

B

are positive constants.

In order to analyze the errors of the difference schemes (41) and (43), by repeatedly reducing the step size from

1 / 40

1 / 200,

the maximum values of errors of

u

and

C

obtained by solving Example 1 are given in tables 1 and 2.

Table 1. The numerical errors using the difference approximation (41) with $α = 0.9,$ $γ = 0.8,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.1,$ $B = - 8,$ $t = 1.$

$h_{x}$	$h_{y}$	$τ$	${∥ e^{k} ∥}_{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 (43) with $α = 0.9,$ $γ = 0.9,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.05,$ $B = - 8,$ $t = 1.$

$h_{x}$	$h_{y}$	$τ$	${∥ ε^{k} ∥}_{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

Let

\begin{array}{rcl} {∥ e^{k} ∥}_{max} = {max}_{1 ⩽ i ⩽ M_{x} - 1, 1 ⩽ j ⩽ M_{y} - 1} | e_{i, j}^{k} |, \end{array}

(67)

\begin{array}{rcl} {∥ ε^{k} ∥}_{max} = {max}_{1 ⩽ i ⩽ M_{x} - 1, 1 ⩽ j ⩽ M_{y} - 1} | ε_{i, j}^{k} | . \end{array}

(68)

From tables 1–2, we can conclude that the errors of the difference scheme (41) and (43) satisfy equations (65), (66).

4.3. Comparison with other method

4.3.1. Example 2

\begin{array}{rcl} \frac{\partial u}{\partial t} = - \frac{\partial u}{\partial x} - \frac{\partial u}{\partial y} + F_{α} \frac{\partial^{α} u}{\partial y^{α}} + f (x, y, t), 0 < α < 1, \end{array}

(69)

with boundary conditions:

\begin{array}{rcl} u (x, y, 0) = 0, 0 ⩽ x ⩽ 1, 0 ⩽ y ⩽ 1, \end{array}

(70)

\begin{array}{rcl} u (0, y, t) = u (1, y, t) = 0, 0 ⩽ y ⩽ 1, 0 ⩽ t ⩽ T, \end{array}

(71)

\begin{array}{rcl} u (x, 0, t) = u (x, 1, t) = 0, 0 ⩽ x ⩽ 1, 0 ⩽ t ⩽ T, \end{array}

(72)

where

\begin{array}{rcl} \begin{array}{l} f (x, y, t) = 2 t x (x {(1 - x)}^{2} + t (1 - 3 x + 2 x^{2})) y^{2} \\ \times {(1 - y)}^{2} + t^{2} x^{2} {(1 - x)}^{2} (2 y (1 - 3 y + 2 y^{2}) \\ - F_{α} q (y)), \end{array} \end{array}

(73)

\begin{array}{rcl} q (y) = \frac{Γ (3) y^{2 - α}}{Γ (3 - α)} - 2 \frac{Γ (4) y^{3 - α}}{Γ (4 - α)} + \frac{Γ (5) y^{4 - α}}{Γ (5 - α)} . \end{array}

(74)

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

α,

parameter

F_{α}

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} = τ = 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} = τ = 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

α

is between 1 and 2, the initial condition

u (x, y, 0) = ϕ (x, y)

is not sufficient, but another initial condition

\frac{\partial u (x, y, 0)}{\partial t} = φ (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.$

$(α, F_{α})$	$h_{x} = h_{y} = τ$	$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

x

and

y

direction, respectively. Figures 6(c) and (d) are the three-dimensional distribution of concentration

C

with respect to time

t

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 $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 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 .$

Full size|PPT slide

Figures 7 and 8 demonstrate the effect of fractional-order constitutive models, i.e. models (4) and (6), the fractional derivatives (

0.8 ⩽ α ⩽ 1.0,

0.7 ⩽ γ ⩽ 1.0

) on the velocity and concentration boundary layers. In figure 7, the increase in fractional-order

α

results in a decrease in horizontal velocity

u .

This is consistent with [27, 28]. The higher the fractional order

α

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

γ .

Overall, the increase in fractional-order

γ

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

γ

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 $α$ for $t = 1,$ $\tilde{R e} = 4,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Full size|PPT slide

Figure 8. Comparison of $C$ with different $γ$ for $t = 1,$ $α = 0.8,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Full size|PPT slide

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 ⩽ \tilde{R e} ⩽ 2.

Tabi et al [58] showed that when

\tilde{R e} = 3,

4,

5,

the flow velocity of fluid increased significantly with the increase of

\tilde{R e} .

In addition, the range of Reynolds numbers were chosen as

\tilde{R e} = 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 ⩽ \tilde{R e} ⩽ 10.

In figure 9(a), as

\tilde{R e}

increases, the flow viscous coefficient

{\tilde{μ}}_{α}

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 $\tilde{R e}$ for $t = 1,$ $α = 0.8,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Full size|PPT slide

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

\tilde{S c}

increases, the molecular diffusion coefficient

{\tilde{D}}_{γ}

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 $\tilde{S c}$ for $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Full size|PPT slide

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 < 0

corresponds to the exudation, and

k > 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,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $n_{1} = 2,$ $n_{2} = 3.$

Full size|PPT slide

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}}_{γ} \frac{\partial^{γ} C}{\partial y^{γ}} .

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}

(m^{2} 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.

Full size|PPT slide

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 $α$ brings about a decrease in the horizontal velocity.
	2. Anomalous diffusion weakens the mass transfer; therefore, the concentration decreases by increasing the fractional derivative $γ .$
	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 .$

References

List( Publishing order | Descend order by publishing year | Descend order by cited within ) Chart analysis

1	Fakhar K Cheng Y Ji X 2006 Lie symmetry analysis and some new exact solutions for rotating flow of a second-order fluid on a porous plate Int. J. Eng. Sci. 44 889 896 https://doi.org/10.1016/j.ijengsci.2006.05.004 Cited in this article [1]

2	Moatimid G M Hassan M A 2012 Viscous potential flow of electrohydrodynamic Kelvin–Helmholtz instability through two porous layers with suction/injection effect Int. J. Eng. Sci. 54 15 26 https://doi.org/10.1016/j.ijengsci.2012.01.002 Cited in this article [1]

3	Maity S 2014 Thermocapillary flow of thin liquid film over a porous stretching sheet in presence of suction/injection Int. J. Heat Mass Transfer 70 819 826 https://doi.org/10.1016/j.ijheatmasstransfer.2013.11.070 Cited in this article [1]

4	Aziz A Aziz T 2012 MHD flow of a third grade fluid in a porous half space with plate suction or injection: an analytical approach Appl. Math. Comput. 218 10443 10453 https://doi.org/10.1016/j.amc.2012.04.006 Cited in this article [1]

5	Idowu A S Falodun B O 2020 Variable thermal conductivity and viscosity effects on non-Newtonian fluids flow through a vertical porous plate under Soret–Dufour influence Math. Comput. Simul. 177 358 384 https://doi.org/10.1016/j.matcom.2020.05.001 Cited in this article [1]

6	Attia H A 2006 On the effectiveness of uniform suction and injection on unsteady rotating disk flow in porous medium with heat transfer Comput. Mater. Sci. 38 240 244 https://doi.org/10.1016/j.commatsci.2006.02.006 Cited in this article [1]

7	Sheremet M A Rosca N C Rosca A V Pop I 2018 Mixed convection heat transfer in a square porous cavity filled with a nanofluid with suction/injection effect Comput. Math. Appl. 76 2665 2677 https://doi.org/10.1016/j.camwa.2018.08.069 Cited in this article [1]

8	Sheikholeslami M 2017 Influence of Lorentz forces on nanofluid flow in a porous cylinder considering Darcy model J. Mol. Liq. 225 903 912 https://doi.org/10.1016/j.molliq.2016.11.022 Cited in this article [1]

9	Sheikholeslami M Seyednezhad M 2017 Nanofluid heat transfer in a permeable enclosure in presence of variable magnetic field by means of CVFEM Int. J. Heat Mass Transfer 114 1169 1180 https://doi.org/10.1016/j.ijheatmasstransfer.2017.07.018 Cited in this article [1]

10	Sheikholeslami M Ganji D 2017 Forced convection in existence of Lorentz forces in a porous cavity with hot circular obstacle using nanofluid via Lattice Boltzmann method J. Mol. Liq. 246 103 111 https://doi.org/10.1016/j.molliq.2017.09.053 Cited in this article [1]

11	Ellahi R Shivanian E Abbasbandy S Hayat T 2015 Analysis of some magnetohydrodynamic flows of third order fluid saturating porous space J. Porous Media 18 89 98 https://doi.org/10.1615/JPorMedia.v18.i2.10 Cited in this article [1]

12	Hassan H Rashidi M 2014 An analytic solution of micropolar flow in a porous channel with mass injection using homotopy analysis method Int. J. Numer. Methods Heat Fluid Flow 242 419 437 https://doi.org/10.1108/HFF-08-2011-0158 Cited in this article [1]

13	Eid M Mahny K Dar A 2020 Numerical study for Carreau nanofluid flow over a convectively heated nonlinear stretching surface with chemically reactive species Physica A 540 123063 https://doi.org/10.1016/j.physa.2019.123063 Cited in this article [1]

14	Arasteh H Mashayekhi R Bahiraei M 2019 Optimal arrangements of a heat sink partially filled with multilayered porous media employing hybrid nanofluid J. Therm. Anal. Calorim. 137 1045 1058 https://doi.org/10.1007/s10973-019-08007-z Cited in this article [1]

15	Saikrishnan P Roy S 2003 Non-uniform slot injection (suction) into water boundary layers over (i) a cylinder and (ii) a sphere Int. J. Eng. Sci. 41 1351 1365 https://doi.org/10.1016/S0020-7225(03)00043-0 Cited in this article [1]

16	Roy S Saikrishnan P 2003 Non-uniform slot injection (suction) into steady laminar water boundary layer flow over a rotating sphere Int. J. Heat Mass Transfer 46 3389 3396 https://doi.org/10.1016/S0017-9310(03)00137-6 Cited in this article [1]

17	Roy S Saikrishnan P 2004 Non-uniform slot injection (suction) into water boundary layer flow past yawed cylinder Int. J. Eng. Sci. 42 2147 2157 https://doi.org/10.1016/j.ijengsci.2003.12.008 Cited in this article [1]

Kumari

Nath

2009 Natural convection from a vertical cone in a porous medium due to the combined effects of heat and mass diffusion with non-uniform wall temperature/concentration or heat/mass flux and suction/injection Int. J. Heat Mass Transfer 52 3064 3069

https://doi.org/10.1016/j.ijheatmasstransfer.2008.10.037

Cited in this article [1]

Ganapathirao

Ravindran

Pop

2013 Non-uniform slot suction (injection) on an unsteady mixed convection flow over a wedge with chemical reaction and heat generation or absorption Int. J. Heat Mass Transfer 67 1054 1061

https://doi.org/10.1016/j.ijheatmasstransfer.2013.08.016

Cited in this article [1]

20	Ravindran R Ganapathirao M 2013 Non-uniform slot suction/injection into mixed convection boundary layer flow over vertical cone Appl. Math. Mech. (Engl. Ed.) 34 1327 1338 https://doi.org/10.1007/s10483-013-1748-7 Cited in this article [1]

21	Li B Yang Y Chen X 2019 A power-law liquid food flowing through an uneven channel with non-uniform suction/injection Int. J. Heat Mass Transfer 144 118639 https://doi.org/10.1016/j.ijheatmasstransfer.2019.118639 Cited in this article [1]

22	Bai Y Huo L Zhang Y 2019 Flow, heat and mass transfer of three-dimensional fractional Maxwell fluid over a bidirectional stretching plate with fractional Fourier's law and fractional Fick's law Comput. Math. Appl. 78 2831 2846 https://doi.org/10.1016/j.camwa.2019.04.027 Cited in this article [3]

23	Mondal S Cassano A Conidi C 2016 Modeling of gel layer transport during ultrafiltration of fruit juice with non-Newtonian fluid rheology Food Bioprod. Process. 100 72 84 https://doi.org/10.1016/j.fbp.2016.06.012 Cited in this article [3]

24	Sheikh N A Ali F Khan I 2018 A modern approach of Caputo-Fabrizio time-fractional derivative to MHD free convection flow of generalized second-grade fluid in a porous medium Neural Comput. Appl. 30 1865 1875 https://doi.org/10.1007/s00521-016-2815-5 Cited in this article [1]

25	Carrera Y Avila-de La Rosa G Vernon-Carter E J 2017 A fractional-order Maxwell model for non-Newtonian fluids Physica A 482 276 285 https://doi.org/10.1016/j.physa.2017.04.085 Cited in this article [1]

26	Sun H Zhang Y Wei S Zhu J Chen W 2018 A space fractional constitutive equation model for non-Newtonian fluid flow Commun. Nonlinear Sci. Numer. Simul. 62 409 417 https://doi.org/10.1016/j.cnsns.2018.02.007 Cited in this article [2]

27	Pan M Zheng L Liu F 2016 Lie group analysis and similarity solution for fractional Blasius flow Commun. Nonlinear Sci. Numer. Simul. 37 90 101 https://doi.org/10.1016/j.cnsns.2016.01.010 Cited in this article [4]

28	Li B Liu F 2020 Boundary layer flows of viscoelastic fluids over a non-uniform permeable surface Comput. Math. Appl. 79 2376 2387 https://doi.org/10.1016/j.camwa.2019.11.003 Cited in this article [2]

29	Xu Z Chen W 2013 A fractional-order model on new experiments of linear viscoelastic creep of Hami Melon Comput. Math. Appl. 66 677 681 https://doi.org/10.1016/j.camwa.2013.01.033 Cited in this article [1]

30	Simpson R Jaques A Nuñez H 2013 Fractional calculus as a mathematical tool to improve the modeling of mass transfer phenomena in food processing Food Eng. Rev. 5 45 55 https://doi.org/10.1007/s12393-012-9059-7 Cited in this article [2]

31	Mehdi D Mostafa A 2018 A finite difference/finite element technique with error estimate for space fractional tempered diffusion-wave equation Comput. Math. Appl. 75 2903 2914 https://doi.org/10.1016/j.camwa.2018.01.020 Cited in this article [1]

32	Liu N Liu Y Li H 2018 Time second-order finite difference/finite element algorithm for nonlinear time-fractional diffusion problem with fourth-order derivative term Comput. Math. Appl. 75 3521 3536 https://doi.org/10.1016/j.camwa.2018.02.014 Cited in this article [1]

33	Yang X Jiang X 2019 Numerical algorithm for two dimensional fractional Stokes' first problem for a heated generalized second grade fluid with smooth and non-smooth solution Comput. Math. Appl. 78 1562 1571 https://doi.org/10.1016/j.camwa.2019.03.029 Cited in this article [1]

34	Ding H Li C 2018 High-order numerical approximation formulas for Riemann–Liouville (Riesz) tempered fractional derivatives: construction and application (II) Appl. Math. Lett. 86 208 214 https://doi.org/10.1016/j.aml.2018.06.037 Cited in this article [1]

35	Yang Y Li M Shu S 2006 High order schemes based on upwind schemes with modified coefficients J. Comput. Appl. Math. 195 242 251 https://doi.org/10.1016/j.cam.2005.04.071 Cited in this article [1]

36	Feng L Liu F Turner I 2019 Finite difference/finite element method for a novel 2D multi-term time-fractional mixed sub-diffusion and diffusion-wave equation on convex domains Commun. Nonlinear Sci. Numer. Simul. 70 354 371 https://doi.org/10.1016/j.cnsns.2018.10.016 Cited in this article [1]

37	Liu L Feng L Zheng L Liu F 2020 Flow and heat transfer of generalized Maxwell fluid over a moving plate with distributed order time fractional constitutive models Int. Commun. Heat Mass Transfer 116 104679 https://doi.org/10.1016/j.icheatmasstransfer.2020.104679 Cited in this article [1]

38	Hao Z Fan K Cao W Sun Z 2016 A finite difference scheme for semilinear space-fractional diffusion equations with time delay Appl. Math. Comput. 275 238 254 https://doi.org/10.1016/j.amc.2015.11.071 Cited in this article [1]

39	Sayevand K Machado J Moradi V 2019 A new non-standard finite difference method for analyzing the fractional Navier–Stokes equations Comput. Math. Appl. 78 1681 1694 https://doi.org/10.1016/j.camwa.2018.12.016 Cited in this article [1]

40	Liu F Zhuang P Anh V 2007 Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation Appl. Math. Comput. 191 12 20 https://doi.org/10.1016/j.amc.2006.08.162 Cited in this article [1]

41	Zhang M Shen M Liu F 2019 A new time and spatial fractional heat conduction model for Maxwell nanofluid in porous medium Comput. Math. Appl. 78 1621 1636 https://doi.org/10.1016/j.camwa.2019.01.006 Cited in this article [1]

42	Wang X Xu H Qi H 2020 Numerical analysis for rotating electro-osmotic flow of fractional Maxwell fluids Appl. Math. Lett. 103 106179 https://doi.org/10.1016/j.aml.2019.106179 Cited in this article [1]

43	Lorenzo G Zaritzky N E Califano A N 2011 Viscoelastic characterization of fluid and gel like food emulsions stabilized with hydrocolloids Procedia Food Sci. 1 281 286 https://doi.org/10.1016/j.profoo.2011.09.044 Cited in this article [1]

44	Yang X Chen W 2017 A fractional model for time-variant Non-Newtonian flow Therm. Sci. 21 61 68 https://doi.org/10.2298/TSCI160426245Y Cited in this article [1]

45	Siginer D 2014 Stability of non-linear constitutive formulations for viscoelastic fluids SpringerBriefs in Applied Sciences and Technology Cham Springer Cited in this article [1]

46	Qi H Xu M 2007 Unsteady flow of viscoelastic fluid with fractiona Maxwell model in a channel Mech. Res. Commun. 34 210 212 https://doi.org/10.1016/j.mechrescom.2006.09.003 Cited in this article [1]

47	Yang D Zhu K 2010 Start-up flow of a viscoelastic fluid in a pipe with a fractional Maxwell's model Comput. Math. Appl. 60 2231 2238 https://doi.org/10.1016/j.camwa.2010.08.013 Cited in this article [1]

48	Ding X Zhang F Zhang G 2020 Modeling of hydraulic fracturing in viscoelastic formations with the fractional Maxwell model Comput. Geotech. 126 103723 https://doi.org/10.1016/j.compgeo.2020.103723 Cited in this article [1]

49	Podlubny I 2000 Matrix approach to discrete fractional calculus Fractional Calculus Appl. Anal. 3 359 386 Cited in this article [2]

50	Welti-Chanes J Vergara-Balderas F Bermúdez-Aguirre D 2005 Transport phenomena in food engineering: basic concepts and advances J. Food Eng. 67 113 128 https://doi.org/10.1016/j.jfoodeng.2004.05.053 Cited in this article [1]

51	John S Marc L 2002 Osmotic dehydration of foods: mass transfer and modeling aspects Food Rev. Int. 18 305 335 https://doi.org/10.1081/FRI-120016208 Cited in this article [1]

52	Simpson R Ramírez C Birchmeier V Almonacid A 2015 Diffusion mechanisms during the osmotic dehydration of Granny Smith apples subjected to a moderate electric field J. Food Eng. 166 204 211 https://doi.org/10.1016/j.jfoodeng.2015.05.027 Cited in this article [1]

53	Ramírez C Astorga V Nunez H 2017 Anomalous diffusion based on fractional calculus approach applied to drying analysis of apple slices: the effects of relative humidity and temperature J. Food Process. Eng. 40 12549 https://doi.org/10.1111/jfpe.12549 Cited in this article [1]

54	Watanabe H Yahata Y Fukuoka M 2007 The thermodynamic basis for the relative water demand model that describes non-Fickian water diffusion in starchy foods J. Food Eng. 83 130 135 https://doi.org/10.1016/j.jfoodeng.2007.02.009 Cited in this article [1]

55	Núñez H Cavada G Ramírez C 2020 Effect of a moderate electric field on the salting of Atlantic Salmon (Salmo salar): an experimental study and phenomenological understanding Food Res. Int. 137 109475 https://doi.org/10.1016/j.foodres.2020.109475 Cited in this article [1]

56	Liu F Anh V Turner I 2004 Numerical solution of the space fractional Fokker-Planck equation J. Comput. Appl. Math. 166 209 219 https://doi.org/10.1016/j.cam.2003.09.028 Cited in this article [1]

57	Sun Z Wu X 2006 A fully discrete difference scheme for a diffusion-wave system Appl. Numer. Math. 56 193 209 https://doi.org/10.1016/j.apnum.2005.03.003 Cited in this article [1]

58	Tabi C B Ndjawa P Motsumi T G 2020 Magnetic field effect on a fractionalized blood flow model in the presence of magnetic particles and thermal radiations Chaos Solitons Fractals 131 109540 https://doi.org/10.1016/j.chaos.2019.109540 Cited in this article [1]

59	Maiti S Shaw S Shit G C 2020 Caputo-Fabrizio fractional order model on MHD blood flow with heat and mass transfer through a porous vessel in the presence of thermal radiation Physica A 540 123149 https://doi.org/10.1016/j.physa.2019.123149 Cited in this article [1]

60	Sebti I Blanc D Carnet-Ripoche A Saurel R Coma V 2004 Experimental study and modeling of nisin diffusion in agarose gels J. Food Eng. 63 185 190 https://doi.org/10.1016/S0260-8774(03)00299-1 Cited in this article [3]

Acknowledgments

The work is supported by the Fundamental Research Funds for the Central Universities and the National Natural Science Foundation of China (12072024).

RIGHTS & PERMISSIONS

PDF(1474 KB)

256

Accesses

Citation

Detail

Sections

Recommended

Abstract
Key words
Cite this article
Nomenclature
1. Introduction
2. Mathematical formulation
Figure 1. Schematic of viscoelastic incompressible liquid food flowing through a permeable plate.
3. Implicit numerical method
4. Algorithm analysis
4.1. Mesh independence
Figure 2. Comparison of $u$ and $C$ with different space step size $h_{x}$ with $x = 0.8,$ $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{y} = 1 / 50,$ $τ = 1 / 50.$
Figure 3. Comparison of $u$ and $C$ with different space step size $h_{y}$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $τ = 1 / 40.$
Figure 4. Comparison of $u$ and $C$ with different time step size $τ$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $h_{y} = 1 / 50.$
4.2. Accuracy and convergence
4.2.1. Example 1
Figure 5. Numerical solution and exact solution of $u,$ $v$ and $C$ with $h_{x} = 1 / 40,$ $h_{y} = 1 / 50$ and $τ = 1 / 40.$
Table 1. The numerical errors using the difference approximation (41) with $α = 0.9,$ $γ = 0.8,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.1,$ $B = - 8,$ $t = 1.$
Table 2. The numerical errors using the difference approximation (43) with $α = 0.9,$ $γ = 0.9,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.05,$ $B = - 8,$ $t = 1.$
4.3. Comparison with other method
4.3.1. Example 2
Table 3. Comparisons between G-L algorithm and L1 algorithm with $t = 1.$
5. Results and discussion
Figure 6. Comparison of $u$ and $C$ with different $t$ for $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 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 .$
Figure 7. Comparison of $u$ with varying $α$ for $t = 1,$ $\tilde{R e} = 4,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$
Figure 8. Comparison of $C$ with different $γ$ for $t = 1,$ $α = 0.8,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$
Figure 9. Comparison of $u$ with varying $\tilde{R e}$ for $t = 1,$ $α = 0.8,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$
Figure 10. Comparison of $C$ with different $\tilde{S c}$ for $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$
Figure 11. Comparison of $C$ with different $k$ for $y = 0.06,$ $t = 3,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $n_{1} = 2,$ $n_{2} = 3.$
Figure 12. Comparison of experimental data from [60] with our results to test the method and the model.
6. Conclusions
References
Acknowledgments
RIGHTS & PERMISSIONS

Received	Revised	Accepted	Published
2020-11-05	2020-12-27	2021-01-25	2021-04-01
Issue Date
2021-04-12

Please choose a citation manager

Content to export

Abstract

Key words

Cite this article

Nomenclature

1. Introduction

2. Mathematical formulation

Figure 1. Schematic of viscoelastic incompressible liquid food flowing through a permeable plate.

3. Implicit numerical method

4. Algorithm analysis

4.1. Mesh independence

Figure 2. Comparison of $u$ and $C$ with different space step size $h_{x}$ with $x = 0.8,$ $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{y} = 1 / 50,$ $τ = 1 / 50.$

Figure 3. Comparison of $u$ and $C$ with different space step size $h_{y}$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $τ = 1 / 40.$

Figure 4. Comparison of $u$ and $C$ with different time step size $τ$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $h_{y} = 1 / 50.$

4.2. Accuracy and convergence

4.2.1. Example 1

Figure 5. Numerical solution and exact solution of $u,$ $v$ and $C$ with $h_{x} = 1 / 40,$ $h_{y} = 1 / 50$ and $τ = 1 / 40.$

Table 1. The numerical errors using the difference approximation (41) with $α = 0.9,$ $γ = 0.8,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.1,$ $B = - 8,$ $t = 1.$

Table 2. The numerical errors using the difference approximation (43) with $α = 0.9,$ $γ = 0.9,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.05,$ $B = - 8,$ $t = 1.$

4.3. Comparison with other method

4.3.1. Example 2

Table 3. Comparisons between G-L algorithm and L1 algorithm with $t = 1.$

5. Results and discussion

Figure 6. Comparison of $u$ and $C$ with different $t$ for $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 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 .$

Figure 7. Comparison of $u$ with varying $α$ for $t = 1,$ $\tilde{R e} = 4,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 8. Comparison of $C$ with different $γ$ for $t = 1,$ $α = 0.8,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 9. Comparison of $u$ with varying $\tilde{R e}$ for $t = 1,$ $α = 0.8,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 10. Comparison of $C$ with different $\tilde{S c}$ for $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 11. Comparison of $C$ with different $k$ for $y = 0.06,$ $t = 3,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 12. Comparison of experimental data from [60] with our results to test the method and the model.

6. Conclusions

{{custom_sec.title}}

{{custom_sec.title}}

References

{{custom_fnGroup.title_en}}

Footnotes

Acknowledgments

{{custom_ack.title_en}}

RIGHTS & PERMISSIONS

模态框（Modal）标题

Please choose a citation manager

Content to export

Abstract

Key words

Cite this article

Nomenclature

1. Introduction

2. Mathematical formulation

Figure 1. Schematic of viscoelastic incompressible liquid food flowing through a permeable plate.

3. Implicit numerical method

4. Algorithm analysis

4.1. Mesh independence

Figure 2. Comparison of u and C with different space step size hx with x=0.8, t=1, α=0.8, γ=0.7, Re~=4, Sc~=4, Cg=6.18, k=−5×10−3, n1=2, n2=3, hy=1/50, τ=1/50.

Figure 3. Comparison of u and C with different space step size hy with t=1, α=0.8, γ=0.7, Re∼=4, Sc∼=4, Cg=6.18, k=−5×10−3, n1=2, n2=3, hx=1/40, τ=1/40.

Figure 4. Comparison of u and C with different time step size τ with t=1, α=0.8, γ=0.7, Re∼=4, Sc∼=4, Cg=6.18, k=−5×10−3, n1=2, n2=3, hx=1/40, hy=1/50.

4.2. Accuracy and convergence

4.2.1. Example 1

Figure 5. Numerical solution and exact solution of u, v and C with hx=1/40, hy=1/50 and τ=1/40.

Table 1. The numerical errors using the difference approximation (41) with α=0.9, γ=0.8, Fα=0.2, Lγ=0.04, s=1, K=0.1, B=−8, t=1.

Table 2. The numerical errors using the difference approximation (43) with α=0.9, γ=0.9, Fα=0.2, Lγ=0.04, s=1, K=0.05, B=−8, t=1.

4.3. Comparison with other method

4.3.1. Example 2

Table 3. Comparisons between G-L algorithm and L1 algorithm with t=1.

5. Results and discussion

Figure 6. Comparison of u and C with different t for α=0.8, γ=0.7, Re∼=4, Sc∼=4, Cg=6.18, k=−5×10−5, n1=2, n2=3: (a) y=0.2; (b) x=0.4; (c) y=0.2; (d) x=0.4.

Figure 7. Comparison of u with varying α for t=1, Re~=4, k=−5×10−4, n1=2, n2=3.

Figure 8. Comparison of C with different γ for t=1, α=0.8, Re∼=4, Sc∼=4, Cg=6.18, k=−5×10−4, n1=2, n2=3.

Figure 9. Comparison of u with varying Re∼ for t=1, α=0.8, k=−5×10−4, n1=2, n2=3.

Figure 10. Comparison of C with different Sc∼ for t=1, α=0.8, γ=0.7, Re∼=4, Cg=6.18, k=−5×10−4, n1=2, n2=3.

Figure 11. Comparison of C with different k for y=0.06, t=3, α=0.8, γ=0.7, Re∼=4, n1=2, n2=3.

Figure 12. Comparison of experimental data from [60] with our results to test the method and the model.

6. Conclusions

{{custom_sec.title}}

{{custom_sec.title}}

References

{{custom_fnGroup.title_en}}

Footnotes

Acknowledgments

{{custom_ack.title_en}}

RIGHTS & PERMISSIONS

Figure 2. Comparison of $u$ and $C$ with different space step size $h_{x}$ with $x = 0.8,$ $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{y} = 1 / 50,$ $τ = 1 / 50.$

Figure 3. Comparison of $u$ and $C$ with different space step size $h_{y}$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $τ = 1 / 40.$

Figure 4. Comparison of $u$ and $C$ with different time step size $τ$ with $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 3},$ $n_{1} = 2,$ $n_{2} = 3,$ $h_{x} = 1 / 40,$ $h_{y} = 1 / 50.$

Figure 5. Numerical solution and exact solution of $u,$ $v$ and $C$ with $h_{x} = 1 / 40,$ $h_{y} = 1 / 50$ and $τ = 1 / 40.$

Table 1. The numerical errors using the difference approximation (41) with $α = 0.9,$ $γ = 0.8,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.1,$ $B = - 8,$ $t = 1.$

Table 2. The numerical errors using the difference approximation (43) with $α = 0.9,$ $γ = 0.9,$ $F_{α} = 0.2,$ $L_{γ} = 0.04,$ $s = 1,$ $K = 0.05,$ $B = - 8,$ $t = 1.$

Table 3. Comparisons between G-L algorithm and L1 algorithm with $t = 1.$

Figure 6. Comparison of $u$ and $C$ with different $t$ for $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $\tilde{S c} = 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 .$

Figure 7. Comparison of $u$ with varying $α$ for $t = 1,$ $\tilde{R e} = 4,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 8. Comparison of $C$ with different $γ$ for $t = 1,$ $α = 0.8,$ $\tilde{R e} = 4,$ $\tilde{S c} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 9. Comparison of $u$ with varying $\tilde{R e}$ for $t = 1,$ $α = 0.8,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 10. Comparison of $C$ with different $\tilde{S c}$ for $t = 1,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $C_{g} = 6.18,$ $k = - 5 \times 10^{- 4},$ $n_{1} = 2,$ $n_{2} = 3.$

Figure 11. Comparison of $C$ with different $k$ for $y = 0.06,$ $t = 3,$ $α = 0.8,$ $γ = 0.7,$ $\tilde{R e} = 4,$ $n_{1} = 2,$ $n_{2} = 3.$