Applications of Orthogonal Polynomials in Simulations of Mass Transfer Diffusion Equation Arising in Food Engineering

Abstract

In this paper, Chebyshev polynomials—which are ultraspherical in the first and second kind and hence symmetric, while the third and fourth order are not ultraspherical and are hence non-symmetric—are used for the simulation of two-dimensional mass transfer equation arising during the convective air drying processes of food products subject to Robin and Neumann boundary conditions. These simulations are used to improve the quality of dried food products and for prediction of the moisture distributions. The equation is discretized in both temporal and special variables by using the second order finite difference scheme and spectral method based on Chebyshev polynomial with the help of fast Fourier transform on tensor product grid, respectively. A system of algebraic equations is obtained after applying the proposed numerical scheme, which is then solved by an appropriate iterative method. The error analysis of the proposed scheme is provided. Some numerical examples are presented to confirm the numerical efficiency and theoretical justification of the proposed scheme. Our numerical scheme has an exponential rate of convergence, which means that one can achieve a very accurate solution using a few collocation points, as opposed to the other available techniques which are very slow in terms of convergence and consume a lot of time. In order to further validate the accuracy of our numerical method, a comparison is made with the exact solution using different norms.

1. Introduction

Orthogonal polynomials and their symmetry and asymmetry properties play a key role in solving differential equations arising in mathematical modeling of some real world problems. They provide a very natural way to interpret, expand, and solve some very useful differential equations. Classical orthogonal polynomials such as Hermit, Legendre, Laguerre as well as discrete ones, including Krawtchouk, Chebyshev, and many others, have widespread applications in many branched of science and engineering. In mathematical analysis and its applications, these polynomials and their symmetric and non-symmetric forms play a crucial role. They are used as basis functions to expand more complicated and complex functions. Chebyshev polynomials are well known for their application in the solution of many linear, nonlinear, delay, integral, and fractional ordinary and partial differential equations subject to some initial and boundary conditions [1,2,3,4,5]. There are six kinds of Chebyshev polynomial, out of which the third and the fourth kind are non-symmetric, while the first, second, fifth, and sixth order are symmetric. Moreover, the application of Legendre polynomials is also used in solving stochastic differential equations [6,7,8]. In addition to the useful properties of orthogonal polynomials, the axis-symmetric geometry also plays an important role in solving complex problems as it can change three-dimensional models to two-dimensional models, which are comparatively easy to solve.

In numerical simulations, interpolation is a common methodology used in engineering sciences for the estimation of values between the data which is known. If the data is limited, interpolation can still produce some smooth surfaces depending on the choice of interpolating polynomials. For the approximation of complex functions or if the working data set is very large, then traditional polynomial interpolation can produce high rates of error. For this reason, one must be very careful when choosing the appropriate interpolation method. Chebyshev polynomials are the best choice for interpolating large data or approximating complex functions as they are more efficient and can produce some very accurate results while interpolating data, together with fast Fourier transform. This quality of Chebyshev polynomials makes them the first choice for approximation while working with some mathematical models arising in food engineering.

Drying technology deals with liquid removal to produce a solid product. It also helps to modify the dried product during the drying process and to improve the final product to meet the quality criteria. To preserve the product quality with minimum energy cost and environmental impact is one of the main goals of drying process. Developing a new product that needs drying or designing a new dryer is very costly and requires a lot of experimental work. Mathematical models play a key role to reduce both cost and experimental work. To model this complex drying process, that is, the removal of moisture from materials that are wet, requires transport equations for heat and mass transfer using conservation law. These equations are highly nonlinear in nature and involve a lot of parameters that affect the processes of drying, some of which may be hard to quantify or measure.

Diffusion equations, a class of partial differential equations (PDEs), are mostly used to model some physical phenomena, especially those arising in food engineering during drying processes. To extend the shelf life of food products, drying technology is used to reduce the cost of packaging and storing [9]. During the process of drying, the mechanism of mass and heat transfer is involved. Fick’s second law of diffusion is used for the mathematical formulation of such phenomena [10,11]. Fick’s second law in general form is given by:

$$\frac{\partial M}{\partial t}=D\frac{{\partial }^{2}M}{\partial x^2},$$

(1)

where M is the concentration, which depends on x and t. D denotes the diffusion coefficients. Various analytical and numerical methods are used to solve the model Equation (1). For the exact solution, the method involves Laplace transform, Green’s function, and Duhamel’s principle [12,13]. The analytical schemes are normally used for comparison with approximate and numerical methods. The main drawback of analytical methods is that they are dealing with infinite series which need a heavy amount of computational time and therefore are not realistic to be used for a such a real world problem [14]. To reduce the cost of computation time, researchers are looking for numerical approaches which are more reliable and robust. To this end, finite difference and spectral methods are the best choices as they both involve derivatives of operational matrices and are hence very popular among the researchers for solving such types of models [15]. A numerical scheme based on finite difference method for two-dimensional heat and moisture transfer during the drying process to predict the temperature and moisture distributions inside the object subjected to convective boundary conditions and the numerical method for coupled one-dimensional heat and mass transport equations with an effective moisture diffusion coefficient are used in [16,17], respectively.

To increase the performance of the drying process, a finite element procedure is used for the transport phenomena’s mathematical model based on theoretical assumptions occurring under the influences of some important variables of air drying as investigated in [18], and for a cooking of rice and water uptake using a Fickian diffusion equation with axisymmetric conditions as used in [19]. The drying of a thin layer of parchment coffee was conducted under controlled temperatures and for a heat and mass transfer during convective drying of porous bodies was studied in [20,21]. An optimal simulation procedure used complex method for dehydration of potatoes (white) to minimize browning in [22]. A Crank–Nicolson scheme to simulate the mathematical expression of arising moisture and heat transfer was used by [23], while an exact solution for the mass and heat transfer equation was obtained in [24]. A higher-order numerical scheme based on Galerkin continuous method for the approximate solution of reaction-diffusion with advection term in multi dimension was used in [25]. Yang et al. used a numerical technique for two-dimensional mass transfer equation based on spectral method [26], while the mathematical model based on Higgs Boson model using the Legendre Galerkin technique was investigated in [27].

Most recently, to enhance the conductivity of high viscosity and lower concentration, a food gum solution was reported in [28]. Investigation during osmodehydration of cubes of apple based on mathematical assumptions was carried out for the mass transfer equation using the data from experiments, and fitted into three types of mathematical models in [29]. A three-dimensional diffusion model for the osmotic dehydration of melon cube pieces subject to mass flows resistance on surfaces was studied in [30]. A fully coupled solid–fluid approach was used for crust freezing in processed meat in [31], while an artificial network technique for convective drying quince slices of cylindrical shape arising in evaluation of moisture contents during convective drying has been studied in [32]. A simulation technique for drying of food processes based on fluid dynamics was studied in detail in [33]. A comprehensive overview to increase the cook–chill shelf-life using innovative packaging methods was carried out in [34]. For more information about the simulation and evaluation of canned foods and their thermal treatment, we refer the reader to [35]. Some other orthogonal polynomials and their properties have been used by many authors, e.g., for the efficient and robust solution of stochastic differential equation with delay with detail convergence, analysis of different models of epidemic diseases and reaction–diffusion model of pattern formation has been conducted in [36,37,38] in [39], respectively. The application of umbral methods to identify 3-variable Hermite and Hermite–Bernoulli polynomials based on Laguerre polynomial has been investigated in [40,41,42].

Spectral methods are the methods which are famous for their exponential rate of convergence. Those that can achieve the maximum possible numerical approximations while needing less computational time are the most commonly used methods for solving PDEs in recent years. They are global in nature, as compared to finite difference and finite elements which are local in nature and require greater computational costs. In general, for time-dependent PDEs, spectral Galerkin or collocation schemes are used with a low order Runge–Kutta method for the resultant system of ordinary differential equations (ODEs) after discretization. As a result, one can get an unrealistic scheme for solving PDEs by using a higher order method in space and a low order method in time. For smooth PDEs, one can use spectral methods in both space and temporal variables. The disadvantage of spectral methods is that they are mostly used in both space and temporal variables in case of one dimension. This deficiency can be covered by their exponential rate of convergence as one needs a few collocation points to achieve a highly accurate solution. To approximate the mass transfer equation, numerical schemes suffer from multiscale effects in the domain of both time and space. Additionally, during the drying processes of foods, the behavior of the ingredients is not linear and needs some very special attention. In order to reduce the effects of these factors, one needs to use a very accurate and efficient numerical scheme for the approximate solution of these models. Most of the existing methods for solving problems of applied nature, such as food engineering models, are deficient and time consuming. To this end, we use a Chebyshev spectral method in space with fast Fourier transform, while a second order finite difference scheme is used in time to solve the model Equation (1) numerically.

The rest of the paper is organized as follows: Section 2 consists of mathematical model and methodology followed by error analysis in Section 3. Numerical simulations are presented in Section 4, followed by conclusion in Section 5.

2. Mathematical Model and Methodology

In food processing, mathematical models play a very crucial role, and can be analytical, empirical, or numerical. These models can provide insight into any process which is complex in nature or difficult to understand. Analytical models are the models that have an analytical solution, which can be easily obtained without the help of any computer. These models were very effective due to their quality of providing complete information about any system in hand. As food processing models are very complicated, analytical models are too simple to capture the very complex properties of food materials. Empirical models come from observation and experimentation and it is sometimes difficult to find a model to start with as the data or experiments may be wrong. For this reason, numerical models are the best tools to be used during the processing of food materials. In food industries there are many unit operations. In most of these unit operations, a transfer of mass and heat accrues which may be steady or unsteady. The most common examples of these operations are drying, heating, cooling sterilization and many more. Out of all these units, the unit of heat transfer is the most important one as it related to nutritional and hygienic quality of the food product.

During the processing of foods, in applying the heat transfer knowledge one must take into account that this can be very difficult with raw materials. For example, in some cases biological materials are not uniform and the resulting product shapes are not regular, sometimes changing drastically during the heating process. During heating processes the physical properties of many products change, as some are time dependent and some dependent on shear. The other problem which frequently occurs during processing is that heat transfer combines with momentum or mass transfer. Apart from this, the effect of heat on the food product also depends on the combination of both time and temperature. During the temperature changes, some biochemical changes may occur. Therefore, during the designing of food thermal processes, heat transfer must be taken into account. The quality of the product will be largely effected due to overheating and may lose its nutritional value, resulting in the spoiling of products. Due to these complexities, the resultant mathematical models may be solved analytically by using many assumptions for simplification. Due to this issue, the use of numerical techniques for solving these models is a must.

Spectral methods are the most powerful tools to solve differential equations, especially partial differential and integral equations arising in mathematical modeling of any physical phenomena. They are global compared with finite element and finite difference which use local polynomials as approximation functions. The exponential convergence property makes the spectral method superior to others, provided when the data in the given equation is smooth. For the problem with regular geometries, spectral methods are much faster and more accurate than any other traditional numerical scheme for the approximate solution of differential equations [43]. In spectral methods the solution of the unknown function, e.g., for the time dependent PDE, can be written as a linear combination of a basis function. The choice of the basis function depends upon the geometry of the problem at hand. For periodic problems, Fourier series is the best choice, while for non-periodic problems orthogonal polynomials are the best choice for basis function. Spectral methods are mainly classified in three different categories: Gelerkin, Tau, and pseudospectral or collocation method. They are different from each other in term of the choices for the basis function. In the case of Gelerkin method the basis function is chosen in a manner to satisfy the boundary conditions (BCs), while in case of Tau method the test function need not to satisfy the BCs and an extra equation is introduced to enforced the BCs. In the case of collocation method these trial functions are delta function, at some points called collocations, and deal with grid points directly similar to finite difference [44,45]. In all the above spectral methods, collocation technique is the one which is mostly used by researchers to find the approximate solution of many complex problems, and is also the main method in the current research.

The orthogonal basis function can be written as a linear combination as:

$$f_N\left(x\right)=\sum _N^{k=0}{\widehat{f}}_k{\varphi }_k\left(x\right).$$

(2)

where the mostly used basis functions are the functions of sine and cosine. In Chebyshev collocation method the Chebyshev polynomial used is:

$$x_i=cos\left(\frac{\pi i}{N}\right).$$

(3)

In the semi circle they can be seen as equidistant points set ${\pi }_i=\frac{\pi i}{N}$, where the Chebyshev polynomial can be calculated as $T_k\left(x\right)=cos\left(kco{s}^{-1}x\right)$, considering k to be the positive integer, the calculated polynomials are:

$$T_0\left(x\right)=1,$$

(4)

$$T_1\left(x\right)=2x^2-1,$$

(5)

$$⋮$$

(6)

$$T_{n+1}\left(x\right)=2x{T}_n\left(x\right)-T_{n-1}\left(x\right).$$

(7)

where the infinite expansion of Chebyshev can be written as:

$$u(x,t)=\sum _{k=0}^{\infty }{b}_k\left(t\right)T_k\left(x\right).$$

(8)

By taking $c_k$ as 2 for k=0 and 1 for $k\ge 1$ our $b_k\left(t\right)$ be

$$b_k\left(t\right)=\frac{2}{c_k}sin\left(\frac{k\pi }{2}\right)j_k\left(\pi \right)e^{-\left(\pi 2\right)t},$$

(9)

For the computation of the discrete Fourier transform and its inverse, an algorithm is used that is known as the fast Fourier transform. Fourier analysis converts the signal from frequency to the original domain and original to the frequency domain. By decomposing the sequence of values into components that have different frequencies discrete Fourier transform was obtained [46]. The discrete Fourier transform is beneficial but has slow computation. Fast Fourier transforms factorize the discrete Fourier transforms matrix into the product of sparse factors [47]. Let $z_0,\dots ,z_{N-1}$ be complex numbers. Discrete Fourier transform is formulated as:

$$Z_k=\sum _{n=0}^{N-1}{z}_n{e}^{-i2\pi kn/N}k=0,1,\dots ,N-1$$

(10)

where $e^{-i2\pi /N}$ is root of unity. Computing this definition requires the operation of order $N^2$. Fast Fourier transform to evaluate the same result in the order of $NlogN$ operations [48].

Consider the mass transfer equation in two-dimensions with the assumptions that heat transfer is negligible, drying air has constant property, and moisture inside is caused by diffusion only. Under these assumptions, the mathematical model for air drying mass transfer can be written in the form [49]:

$$\frac{\partial M}{\partial t}=D\left[\frac{{\partial }^{2}M}{\partial x^2}+\frac{{\partial }^{2}M}{\partial y^2}\right],$$

(11)

where $(x,y)\in [a,b]\times [c,d]$ with the initial and boundary conditions as

$$M(x,y,0)=M_0,$$

$$\frac{\partial M(a,y,t)}{\partial x}=\frac{\partial M(x,c,t)}{\partial y}=0$$

$$-D\frac{\partial M(b,y,t)}{\partial x}=h_m(M(b,y,t)-M_{air})$$

(12)

$$-D\frac{\partial M(x,d,t)}{\partial y}=h_m(M(x,d,t)-M_{air})$$

Here D represents the diffusion parameter, $h_m$ is the coefficient of mass transfer, and $M_{air}$ along with $M_0$ are time-dependent constants. It must be noted that the model Equation (11), subject to the given initial and boundary conditions given in Equation (12) can be transformed from constant coefficients to variable coefficients, but for simplicity we here consider it with constant coefficient to solve the model equation efficiently. To this end, we use spectral method together with the fast Fourier transform, which speeds up the calculation drastically. In this method, z is introduced as the complex number on the unit circle, and the argument of z is $\varphi$ [50]. Suppose that x is the real part of z and is represented by:

$$x=\frac{1}{2}\left(z+z^{-1}\right)=cos\varphi$$

(13)

Hence the $nth$ polynomial is represented as:

$$P_n\left(x\right)=\frac{1}{2}\left(z^n+z^{-n}\right)=cosn\varphi .$$

(14)

Which generally turns out to be:

$$P_{n+1}\left(x\right)=\frac{1}{2}\left(z^n+z^{-n}\right)\left(z+z^{-1}\right)-\frac{1}{2}\left(z^{n-1}+z^{1-n}\right),$$

(15)

giving the relation as:

$$P_{n+1}\left(x\right)=2x{P}_n\left(x\right)-P_{n-1}\left(x\right),$$

(16)

These polynomials can be uniquely written as a linear combination as:

$$p\left(x\right)=\sum _{n=0}^N{a}_n{P}_n\left(x\right),$$

(17)

resulting as:

$$p\left(\varphi \right)=\sum _{n=0}^N{a}_{n}cosn\varphi .$$

(18)

Forming the self-reciprocal function as:

$$F\left(z\right)=f\left(\frac{1}{2}\left(z+z^{-1}\right)\right),F\left(\varphi \right)=f\left(cos\varphi \right)$$

(19)

By using the coefficients $a_n$ we get:

$$p\left(\varphi \right)=\frac{1}{2\pi }\sum _{k=-N+1}^N{e}^{ik\varphi }{v}_k^{/}=\sum _{n=0}^N{a}_{n}cosn\varphi .$$

(20)

The main idea is to carry out Chebyshev differentiation by fast Fourier transform with collocations points equally distributed on the space of a unit circle.

To apply the spectral method, we deal with the interpolation points as:

$${\varphi }_j=j\pi /N,x_j=cos{\varphi }_j$$

(21)

The derivative of algebraic polynomial interpolation is:

$$w^{\prime }\left(x\right)=\frac{W^{\prime }\left(\varphi \right)}{dx/d\varphi }=\frac{{\sum }_{n=0}^{N}n{a}_{n}sinn\varphi }{\sqrt{1-x^2}},$$

(22)

hence the second derivative is then calculated as:

$$w^{″}\left(x\right)=\frac{-x_j}{{(1-x_j^2)}^{3/2}}{W}_j+\frac{1}{1-x_j^2}{W}_j^{″},$$

(23)

where j is from 0 to N.

To solve the considered equation numerically, finite difference method of second order is used in time derivative t and Chebyshev spectral differentiation with the help of fast Fourier transform on tensor product grid in x and y.

$$\frac{\partial M}{\partial t}=\frac{M_{i+1}-M_{i-1}}{2\Delta t},$$

(24)

similarly, for spatial derivatives we get:

$$\frac{{\partial }^{2}M}{\partial x^2}=w^{″}\left(x\right),$$

(25)

$$\frac{{\partial }^{2}M}{\partial y^2}=w^{″}\left(y\right),$$

(26)

resulting the equations as:

$$\frac{M_{i+1}-M_{i-1}}{2\Delta t}=w^{″}\left(x\right)+w^{″}\left(y\right),$$

(27)

$$M_{i+1}=M_{i-1}+2\Delta t\left[w^{″}\left(x\right)+w^{″}\left(y\right)\right],$$

(28)

where the initial and boundary conditions were also implemented accordingly.

3. Error Analysis

To prove the error analysis of our method we first introduce some useful lemmas

Lemma 1.

Assume that a $(N+1)$-point Gauss, or Gauss–Radau, or Gauss–Lobatto quadrature formula relative to the weight, which is used to integrate the product $y\varphi$, where $y\in H^m\left(I\right)$ with $I:=(-1,1)$ for some $m\ge 1$ and $\varphi \in {\mathcal{P}}_N$. Then there exists a constant C independent of N such that:

$$\left|u\left(t\right)\varphi \left(t\right)dt-{(u,\varphi )}_N\right|\le C{N}^{-m}{\left|u\right|}_{m,N}\left(I\right){\parallel \varphi \parallel }_{L^2\left(I\right)},$$

(29)

where

$${\left|u\right|}_{m,N}\left(I\right)={\left(\sum _{k=min(m,N)}^m{∥u^{\left(k\right)}∥}_{L^2\left(I\right)}^2\right)}^{1/2}$$

and

$${(u,\varphi )}_N=\sum _{k=0}^N{\omega }_{k}u\left(t_k\right)\varphi \left(t_k\right).$$

Lemma 2.

Assume that $u\in H^m\left(I\right)$ and denote $I_{N}y$ the interpolation polynomial associated with the $(N+1)$-point Gauss, or Gauss–Radau, or by Gauss–Lobatto points ${\left\{t_k\right\}}_{k=0}^N$. Then

$$\parallel u-I_N{u\parallel }_{L^2\left(I\right)}\le C{N}^{-m}{\left|u\right|}_{{}_{m,N}\left(I\right)},$$

$$\parallel u-I_N{u\parallel }_{L^{\infty }\left(I\right)}\le C{N}^{3/4-m}{\left|u\right|}_{{}_{m,N}\left(I\right)}.$$

Proof. 

The above estimate is given on p. 289 of [51]. The following estimate

$$\parallel u-I_N{u\parallel }_{L^2\left(I\right)}\le C{N}^{-m}{\left|u\right|}_{{}_{m,N}\left(I\right)}$$

$$\parallel u-I_N{u\parallel }_{L^{\infty }\left(I\right)}\le C{N}^{3/4-m}{\left|u\right|}_{{}_{m,N}\left(I\right)}$$

is also given in [51]. Using the above estimate and the inequality

$${\parallel v\parallel }_{L^{\infty }(a,b)}\le \sqrt{\frac{1}{b-a}+2}{\parallel v\parallel }_{L^2(a,b)}^{1/2}{\parallel v\parallel }_{H^1(a,b)}^{1/2},\forall v\in H^1(a,b).$$

□

Now we want to show that the error between the exact and approximate solution decreases exponentially, that is:

$$\parallel U^{\left(k\right)}-u^{\left(k\right)}{\parallel }_{L^2\left(I\right)}\to 0,asN\to \infty$$

where U is the polynomial of degree N associated with the spectral approximation Equation (19) and C is a constant independent of N. To prove this claim, using the fact that the first kind Chebyshev polynomials ${\phi }_k\left(t\right)$ are the eigenfunctions of the singular Sturm–Liouville of the form [52,53]:

$${\left(1-t^2\right)}^{1/2}{\phi }_k^{\prime }{\left(t\right)}^{\prime }+\frac{k^2}{{\left(1-t^2\right)}^{1/2}}{\phi }_k\left(t\right)=0,t\in (-1,1)$$

(30)

If ${\phi }_k\left(t\right)$ is normalized, that is ${\phi }_k\left(1\right)=1$, then:

$${\phi }_k\left(t\right)=cos\left(m\theta \right),\theta =co{s}^{-1}t.$$

(31)

Let $P_N\left(t\right)$ denote the $N^{th}$ order Chebyshev polynomials of degree $\le N$. These polynomials are well known for their orthogonality on the interval $[-1,1]$ with respect to the weight function $\omega \left(t\right)$, that is:

$${\int }_{-1}^1{\phi }_k\left(t\right){\phi }_l\left(t\right)\omega \left(t\right)dt=0,k\ne l.$$

(32)

Such system is complete in the space $L_{\omega }^2(-1,1)$, which is space of function $f\left(t\right)$ on $(-1,1)$ of the form:

$${\int }_{-1}^1{\left|f\left(t\right)\right|}^2\omega \left(t\right)dt<\infty .$$

(33)

The space $L_{\omega }^2(-1,1)$ is Hilbert space with respect to the inner product by lemma 1:

$$\langle f,g\rangle ={\int }_{-1}^{1}f\left(t\right)g\left(t\right)\omega \left(t\right)dt,$$

(34)

where the norm:

$${\parallel f\parallel }_{L_{\omega }^2(-1,1)}={\left({\int }_{-1}^1{\left|f\left(t\right)\right|}^2\omega \left(t\right)dt\right)}^{1/2}.$$

(35)

Now, if the solution $u\left(t\right)\in L_{\omega }^2(-1,1)$, then one can write the unknown as series of the function ${\phi }_k\left(t\right)$ of the form:

$$u\left(x\right)=\sum _{k=0}^{\infty }{\phi }_k\left(t\right)\widehat{u_k},$$

(36)

where the coefficient of expansion $\widehat{u_k}$ is given by:

$$\frac{{\int }_{-1}^{1}u\left(x\right){\phi }_k\left(t\right)\omega \left(t\right)dt}{\parallel {\phi }_k{\left(t\right)\parallel }_{L_{\omega }^2(-1,1)}^2}.$$

(37)

The numerical solution $\widehat{u_k}$ can be written as the truncated Cheybeshev series form:

$$p_N\widehat{u_k}=\sum _{k=0}^N\widehat{u_k}{\phi }_k\left(t\right).$$

(38)

where $p_N$ is orthogonal projection of u on $P_N$ in the inner product

$${\langle P_N\widehat{u},\psi \rangle }_{\omega }={\langle y,\psi \rangle }_{\omega },\forall \psi \in P_N.$$

(39)

The completeness of the ${\phi }_k$ is equivalent to the property by Lemma 2:

$$\parallel u-p_N\widehat{u}{\parallel }_{\omega }\to 0,asN\to \infty ,\forall u\in L_{\omega }^2(-1,1).$$

(40)

One can see that Chebyshev differentiation method can be viewed as the hp version of the finite element method. Following this process, we can approach the error estimation of the form $O\left(e^{-cN}\right)$, which is the exponential order of convergence.

4. Numerical Examples

To confirm the accuracy and the exponential rate of convergence of the proposed scheme, two numerical examples are presented in this section. Highly accurate numerical solutions are obtained with the help of the proposed numerical scheme. Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 are obtained for different values of t. A comparison is made with the exact solution, which confirms our theoretical justification for exponential order of convergence. In order to further validate our exponential order of convergence, we provide the error between the exact and approximate solution in

Table 1. Example 1: Error behavior between exact and numerical solution.

Error	$\mathit{t}=0$	$\mathit{t}=0.333$	$\mathit{t}=0.666$	$\mathit{t}=1$
$L_{\infty }$	0	$5.3291e^{-15}$	$1.4211e^{-14}$	$2.1316e^{-14}$
$L_2$	0	$1.16e^{-14}$	$2.8732e^{-14}$	$7.236e^{-14}$
$L_1$	0	$3.7637e^{-14}$	$1.1324e^{-13}$	$3.076e^{-13}$

and

Table 2. Example 2: Error behavior between exact and numerical solution.

Error	$\mathit{t}=0$	$\mathit{t}=0.333$	$\mathit{t}=0.666$	$\mathit{t}=1$
$L_{\infty }$	0	$7.7061e^{-5}$	$3.9541e^{-5}$	$2.0301e^{-5}$
$L_2$	0	$9.8563e^{-4}$	$5.057e^{-4}$	$2.5965e^{-4}$
$L_1$	0	$0.0318$	$0.0071$	$0.0036$

for both example 1 and example 2, respectively.

4.1. Example 1

The two-dimensional mass transfer equation is considered as:

$$\frac{\partial M}{\partial t}=\frac{{\partial }^{2}M}{\partial x^2}+\frac{{\partial }^{2}M}{\partial y^2},$$

(41)

with the initial and boundary conditions as:

$$M(x,y,0)=e^{(x+y)},$$

$$\frac{\partial M(-1,y,t)}{\partial x}=e^{y-1+2t}$$

$$\frac{\partial M(x,-1,t)}{\partial y}=e^{x-1+2t}$$

(42)

$$M(1,y,t)+\frac{\partial M(1,y,t)}{\partial x}=2e^{y+1+2t}$$

$$M(x,1,t)+\frac{\partial M(x,1,t)}{\partial x}=2e^{x+1+2t}$$

Analytical solution is given by:

$$M(x,y,t)=e^{(x+y+2t)},$$

The comparison between the approximate and exact solution is given in Figure 1, Figure 2, Figure 3 and Figure 4.

4.2. Example 2

The example constructed here is also two-dimensional and is presented as:

$$\frac{\partial M}{\partial t}=\frac{{\partial }^{2}M}{\partial x^2}+\frac{{\partial }^{2}M}{\partial y^2},$$

(43)

with the initial and boundary conditions as

$$M(x,y,0)=cos(x+y),$$

$$\frac{\partial M(-1,y,t)}{\partial x}=-e^{-2t}sin(y-1)$$

$$\frac{\partial M(x,-1,t)}{\partial y}=-e^{-2t}sin(x-1)$$

(44)

$$M(1,y,t)+\frac{\partial M(1,y,t)}{\partial x}=e^{-2t}\left(cos(y+1)sin(y+1)\right)$$

$$M(x,1,t)+\frac{\partial M(x,1,t)}{\partial x}=e^{-2t}\left(cos(x+1)sin(x+1)\right)$$

Analytical solution of Equation (43) is

$$M(x,y,t)=e^{(-2t)}cos(x+y),$$

Using the proposed scheme, simulations are made for different time t. A comparison with exact solution is shown in Figure 5, Figure 6, Figure 7 and Figure 8.

5. Conclusions

It has been the aim of this paper to provide a robust, efficient, and accurate approximate solution for the two-dimensional mass transfer equation arising in food engineering. A Chebyshev spectral method along with fast Fourier transform is applied in space coordinates, while the temporal part is discretized using finite difference scheme. It is shown both theoretically as well numerically with the help of some numerical examples that the proposed scheme has exponential order of convergence. To further validate our results, we also make comparison between the exact solution and approximate solution. In addition to this, very few collocation points are needed to get an accurate approximate solution. We find the numerical solution has very good agreement with the exact solution, as our model is based on mathematical formulation of drying of food engineering processes and the geometry of most of the food products is spherical or cylindrical. One of our future tasks is to solve these models in cylindrical and spherical coordinate systems. We also intend to extend the same model into three dimensions.