Legendre Wavelet Modified Petrov–Galerkin Method in Two-Dimensional Moving Boundary Problem

1 Introduction

The applications of moving boundary problems frequently appear in physical, engineering and industrial problems like the melting and solidification of alloys [1], [2], [3], [4], casting and welding of metals and alloys [5], [6], growth of crystals from melts [7], [8], [9], thermal energy storage [10], [11], etc. The solution of moving boundary problems is one of the challenging tasks due to the existence of nonlinear interface condition which is moving with time. Therefore, very few analytical solutions exist only in one-dimensional case with simple initial and boundary conditions. Several numerical methods were used for solving moving boundary problems such as front-tracking methods, front-fixing and fixed domain. These methods are discussed in detail by Crank [12]. The heat balance integral method [13], variable space grid method [14], enthalpy method [15] and coordinate transform method [16] are established for solving one-dimensional single-phase moving boundary problems. Goodman and Shea [17] applied the heat balance integral method to the two-phase problems of melting of a finite slab. Further, to improve the accuracy of heat balance integral method Noble [18] proposed a spatial subdivision scheme in which quadratic profiles were used in each subregion. After that, the Noble’s scheme was modified by Bell [19]. Poots [20] implemented the integral method in two-dimensional moving boundary problem. He considered the two-dimensional temperature profile when liquid is inside a square prism. Rasmussen [21] solved a weak two-dimensional problem choosing a linear profile for the temperature distribution in one of the direction at all times. Gupta and Kumar [22] suggested two methods; in one method a one-dimensional quadratic profile was chosen at preselected values of the other variable and in the second method a piecewise linear profile was assumed. Gupta and Banik [23], [24] proposed the constrained integral method to solve the same problem when the surface was subjected to the boundary condition of the first kind. Instead of choosing two-dimensional temperature profile, he used the temperature profile for constant x. A general implicit source-based enthalpy method was used by Swaminathan and Voller [25] to solve two-dimensional solidification problems with a mushy zone. The multi-dimension moving boundary problems are much harder to solve than the one-dimensional problems.

In last two decades, the application of finite difference method [26], finite element method [27], [28], [29] and boundary element method [30], [31] were used extensively for solving moving boundary problems. On the other hand, several approximate methods and computational techniques were also proposed to handle the nonlinearity of moving boundary problems. For example, Kharab [32] used an approximate method namely, Karman–Polhausan method for the two-dimensional parabolic free boundary problem. Rai and Singh [33] provided a numerical solution of the two-phase moving boundary problem in a finite domain. Yagit [34] proposed an approximate analytical solution and also proposed numerical solutions for a two-dimensional Stefan problem by using a linear perturbation method. Slota [35], [36] employed the genetic algorithm for solving the two-dimensional two-phase inverse Stefan problem. Liu et al. [37] used a quasi-reversible method to identify a moving boundary for two-dimensional heat conduction in a bounded domain without using initial temperature. Chen and Raju [38] developed the coupled finite element and meshless local Petrov–Galerkin method for two-dimensional potential problems. In combined finite element and meshless local Petrov–Galerkin method, the trial functions are taken as the sum of shape function of the each element used in traditional finite element method. The test functions are chosen to be the variations on the sum of shape function as is customary in the conventional finite element method.

The other two most attractive methods were the level set method and phase-field method, and these methods were widely used for solving moving boundary problems. The level set method was firstly, developed by Osher and Sethian [39]. The application of level set method was already generalised to many problems [40], [41], [42]. In this method, the interface position is captured by the level set function and set as its zero level, and it is advected by the introduction of a hyperbolic equation into the governing set of equations. The simple phase-field method was developed by Kobayashi [43] for one-component melt growth including anisotropy. Further, this method is extended by Wheeler et al. [44] for binary alloys. In phase-field method, the time evolution equation is solved by the newly introduced phase-field variable, which equals a fixed constant in each phase. In this approach, the sharp interface is replaced by a diffuse interface with a finite thickness, where the phase-field variable varies rapidly and smoothly. Javierre et al. [45] solved one-dimensional Stefan problems by using the moving grid, level set and phase field methods and applied similarity solutions to verify the numerical results. Takaki et al. [46] used two-dimensional phase-field simulation for directional solidification of Al-Cu alloy. Recently, Yadav et al. [47] solved a one region inward solidification problem whose one surface is kept under most generalised boundary condition and used the finite element Legendre wavelets Galerkin method. Upadhyay et al. [48] applied Legendre wavelet Galerkin method in free boundary problems of melting and solidification, and found that error increases as Stefan number increases. Yadav et al. [49] developed the new approach of the spectral Galerkin method with Legendre wavelets as basis functions for the solution of solidification of the binary eutectic system. Furthermore, they have proved the well-posedness of the model and the stability analysis of the method. The Legendre wavelet with different methods (Galerkin, collocation, spectral) gives better results as compared to other wavelets and the solution obtained by these methods converge to the exact solution for large values of shifting parameters [50], [51], [52].

In this paper, we studied a two-dimensional moving boundary problem arising during melting of solid whose one surface is kept under most generalised boundary conditions and the other two surfaces insulated. We developed the Legendre wavelet modified Petrov–Galerkin method for the solution of two-dimensional moving boundary problem. In the process of the original differential equation, the test function chosen as the linear combination of two-dimensional Legendre wavelet basis functions and trial function are taken as the two-dimensional Legendre wavelet basis functions. For the initial and boundary conditions, we chose the trial function as same and the test function as the first element of the Legendre wavelet basis function. The error minimised by considering that the residual obtained from the differential equation is orthogonal to Legendre wavelet basis function and the residuals obtained from initial and boundary conditions are orthogonal to first Legendre wavelet basis function. Further, we have proved the convergence analysis of this method and also studied the effect of the parameters such as Kirpichev number, Biot number, Predvoditelev number and Fourier number over moving layer thickness.

2 Formulation of Mathematical Model

Let Ω={(x₁, y₁, t)|0≤x₁≤l, 0≤y₁≤l, 0<t<∞} be a domain occupied by solid initially at a temperature g₀(x₁, y₁). The surface y₁=0 is heated by imposing on it the boundary condition of first kind or second kind or third kind. Namely, one might assume either a constant temperature T_s>T_m or a constant heat flux q or a constant heat transfer co-efficient α. The other two surfaces are insulated. At time t>0, the melting starts and liquid/solid interface s(x₁, t) moves in y₁-direction. Therefore, the domain Ω is divided into two regions. The first domain D={(x₁, y₁, t)|0≤x₁≤l, 0≤y₁≤s(x₁, t)} consisting of melt and D′={(x₁, y₁, t)|0≤x₁≤l, s(x₁, t)≤y₁≤l} consisting of solid. In solid region the temperature T_m is constant throughout and the temperature in liquid region is unknown. Hence, the problem is a one-phase problem and the mathematical formulation for moving boundary problem is described by the differential equations:

where T is the temperature and a is the thermal diffusivity. The initial and boundary conditions are

where,

1. A₀=0, B₀=1, f₀(t)=T_s,

2. A₀=K₀, B₀=0, f₀(t)=q,

3. A₀=K₀, B₀=α, f₀(t)=αT_m,

are defined in the first, second and the third kind boundary conditions, respectively.

The liquid/solid interface conditions are

3 Legendre Wavelets Modified Petrov–Galerkin Method

Let us assume that the approximation of θ(x, y, Fo) is a multiplication of e^−bFo and the linear combination of two-dimensional Legendre wavelets basis functions, i.e.

where C and ψ(x, y) are the column vector of order 2^k−12^k′−1MM′×1 which is given by

C=[c1010,⋯,c101M′−1, c1020⋯c102M′−1, c102k′−10⋯c102k′−1M′−1, ⋯,c1M−110, ⋯, c1M−11M′−1, c1M−120, ⋯, c1M−12M′−1, ⋯, c1M−12k′−10, ⋯,c1M−12k′−1M−1, c2010, ⋯, c201M′−1, c2020, ⋯, c202M′−1, ⋯, c202k′−10, ⋯,c202k′−1M′−1, c2M−110, ⋯, c2M−11M′−1, c2M−120, ⋯, c2M−12M′−1, ⋯,c2M−12k′−10, ⋯, c2M−12k′−1M′−1, ⋯, c2k−1010, ⋯, c2k−101M′−1, c2k−1020, ⋯,c2k−102M′−1, ⋯, c2k−102k′−10, ⋯, c2k−1M−12k′−1M′−1]T,

and

ψ=[ψ1010, ⋯, ψ101M′−1, ψ1020⋯ψ102M′−1, ψ102k′−10⋯ψ102k′−1M′−1, ⋯,ψ1M−110, ⋯, ψ1M−11M′−1, ψ1M−120, ⋯, ψ1M−12M′−1, ⋯,ψ1M−12k′−10, ⋯, ψ1M−12k′−1M−1, ψ2010, ψ201M′−1, ψ2020, ⋯,ψ202M′−1, ⋯, ψ202k′−10, ⋯, ψ202k′−1M′−1, ψ2M−110, ⋯,ψ2M−11M′−1, ψ2M−120, ⋯, ψ2M−12M′−1, ⋯, ψ2M−12k′−10, ⋯,ψ2M−12k′−1M′−1, ⋯, ψ2k−1010, ⋯, ψ2k−101M′−1, ⋯, ψ2k−1020, ⋯,ψ2k−102M′−1, ⋯, ψ2k−102k′−10, ⋯, ψ2k−1M−12k′−1M′−1]T.

The elements of column vector ψ are the two-dimensional Legendre wavelets which is defined in the interval [0, 1]×[0, 1] as follows [53]:

where m=0, 1, 2, ……, M−1, m′=0, 1, 2, ….., M′−1, n=1, 2…., 2^k−1 and n′=1, 2…., 2^k′−1, n^=2n−1n^′ =2n′−1.P_m(x) and P_m′ (y) are well-known Legendre polynomials of order m and m′, respectively, and are given as:

P0(x)=1, P1(x)=x, (m+1)Pm+1(x)=(2m+1)Pm(x)−mPm−1(x),P0(y)=1, P1(y)=y, (m′+1)Pm′+1(y)=(2m′+1)Pm′(y)−m′Pm′−1(y),

respectively.

The basis function e^−bFoψ(x, y) is taken in such a way that the Euclidean norm of basis function is a unit. This implies that b=1. Differentiate (14) with respect to Fo, two times with respect to x and y and using the following properties [54]:

∂ψ(x, y)∂x=Dxψ(x, y),∂ψ(x,y)∂y=Dyψ(x, y),

we get,

where D_x and D_y are the operational matrices of differentiation of order 2^k−12^k′−1MM′×2^k−12^k′−1MM′ which are defined as follows:

Dx=(DN⋯NND⋯N⋮⋮⋱⋮NN⋯D), Dy=(D′N⋯NND′⋯N⋮⋮⋱⋮NN⋯D′),

where N is MM′×MM′ null matrix and D, D′ are MM′×MM′ matrices given by

D=(EO⋯OOE⋯O⋮⋮⋱⋮OO⋯E), D′=(FO′⋯O′O′F⋯O′⋮⋮⋱⋮O′O′⋯F).

In which O and O′ are null matrices of order M×M and M′×M′, respectively. E=(e_r,s ) and F=(f_r,s ) be M×M and M′×M′ matrices, respectively. The elements of E and F are defined as follows:

er,s={2k(2r−1)(2s−1), r=2, 3, ⋯M;s=1, 2, …r−1,   r+s=odd0, otherwise.fr,s={2k′(2r−1)(2s−1), r=2, 3, ⋯M′;s=1, 2, …r−1, r+s=odd0, otherwise.

Substituting (16)–(18) in (8)–(12), we obtained the residual R₀(x, y, Fo, C) of differential (8) as

taking the transpose (19), we get

and the residuals of initial and boundary conditions of (9)–(12) are obtained as follows:

As the solid consists of {(x, y, Fo)|0≤x≤1, 0≤y≤1} and λ(x, Fo) varies on y-axis and lies from 0 to 1. So, for minimising the error, we assume that the residual obtained from differential (20) is orthogonal to wavelet basis functions e^−Foψ(x, y), i.e.

After simplifying (26) and using the property

we get

where, Ō is a zero column vector of order 2^k−12^k′−1MM′×1. Equation (28) represents the 2^k−12^k′−1MM′ system of algebraic equations. Further, we assume that the residuals obtained from initial and boundary conditions are orthogonal to the first basis function e^−Foψ₁₀₁₀(x, y) i.e.

Equations (29)–(33) represent the five algebraic equations in unknown C.

Without loss of generality, we take first 2^k−12^k′−1MM′−5 algebraic equations from (28) and remaining five from (29)–(33) that makes the system of algebraic equations of order 2^k−12^k′−1MM′ in unknown vector C and gives the unique solution provided; the determinant is not vanishing. By using Gauss elimination technique, we obtained the vector C in terms of λ. Substituting the value of C in (14), we thus obtained θ in terms of x, y, λ and Fo.

4 Numerical Computation and Discussion

For analysis of model, the particular function g(x, y) is taken as

g(x, y)=1−ygj(x, 0),  j=1, 2

where,

Case (I):

g1(x, 0)=5−2  cos(πx/2)10,

Case (II):

g2(x, 0)=5−2 cos(πx)10.

More details of particular function g(x, y) are given by Gupta and Banik [24]. The particular cases of the generalised boundary condition are taken as

Case 1:

Case 2:

where Ki is the Kirpichev number defined as Ki=qlK0(Ts−Tm) and Pd is the Predvoditelev number defined as Pd=bl2K0(Ts−Tm),b is certain time-dependent coefficient.

Case 3:

where Bi is the Biot number defined as Bi=αlK0.

Substituting these values of parameters defined in (36)–(38) in (11), we get the boundary condition of a first kind, second kind, and third kind, respectively. For numerical computation of this problem, we have used the two-dimensional Legendre wavelets Petrov–Galerkin method defined in Section 3. We have taken M=3, M′=3, k=1, k′=1 for wavelet basis functions. Therefore, we get nine wavelet basis functions. All the numerical computations were done by Matlab (R2012a) (Designtech Systems Ltd., Tal Mulshi, Pune, India). For the validity of the present method, we have plotted graphs (Figs. 1–8) between residuals (obtained from the original differential equation and its associated boundary conditions) and x-axis. We found the effect of an error on moving layer thickness along x and y coordinates, respectively.

Figure 1:

Residual obtained from differential equation in first kind boundary condition at Fo=5 Case (II).

Figure 2:

Residual obtained from boundary condition defined at y=λ in first kind boundary condition at Fo=5 Case (II).

Figure 3:

Residual obtained from boundary condition defined at y=λ in second kind boundary condition at y=λ=Fo=0.5, Ki=Pd=0.01 Case (II).

Figure 4:

Residual obtained from boundary condition defined at y=λ in second kind boundary condition at y=λ=Fo=0.5, Ki=Pd=0.01 Case (II).

Figure 5:

Residual obtained from boundary condition defined at y=λ in second kind case boundary condition at y=λ=Fo=Pd=0.5 Case (II).

Figure 6:

Residual obtained from boundary condition defined at y=λ in second kind boundary condition at y=λ=Fo=Ki=0.5 Case (II).

Figure 7:

Residual obtained from differential equation in third kind boundary condition at Fo=5 Case (II).

Figure 8:

Residual obtained from boundary condition defined at y=λ in third kind boundary condition at Fo=5 Case (II).

From Figure 1, we observed that the error [residual obtained from differential equation (19)] increases as dimensionless moving layer thickness increase at y=0.01, Fo=5 and errors have alternative behaviour when x is increasing for boundary condition of the first kind. It is apparent from Figure 2, the errors [residual obtained from (24)] have an alternate behaviour when x and y are increasing, for boundary condition of the first kind. For boundary condition of the second type, the error increases as λ increases as shown in Figures 3 and 4, We observed that the error decreases as Fo increases.

In Figures 5 and 6, for the second kind boundary condition, the error remains zero on x for different values of Ki and Pd, respectively. In the case of a boundary condition of the third kind, Figure 7 shows that the error [residual obtained from differential equation (19)] increases as y decreases and x increases. In Figure 8, at y=λ the errors [residual obtained from differential equation (24)] tends to zero for boundary condition of the third kind.

For first, second and third kind boundary conditions, we have computed the moving layer thickness numerically by Runge–Kutta method for different values of x and Fo. Figures 9 and 10 show the graph for moving layer thickness for Cases (I) and (II). In these figures, at x=1, the dimensionless moving layer thickness increases as Fourier number (Fo) increases. By the definition of λ and Fourier number Fo, we conclude that moving layer thickness increases as thermal diffusivity of material and time t increase and length of solid l decreases. That means the rate of moving layer thickness is increasing as the volume of solid is decreasing. From these figures, we also observed that the time required for melting of solid is less in the boundary condition of a third kind as compared to the boundary condition of a second kind and less in boundary condition of the second kind in comparison to the first kind.

Figure 9:

Moving layer thickness λ(Fo) with dimensionless time Fo Case (I).

Figure 10:

Moving layer thickness λ(Fo) with dimensionless time Fo Case (II).

In the second kind boundary condition, we studied the effect of parameters such as Predvoditelev number Pd and Kirpichev number Ki on moving layer thickness λ(Fo). In Figures 11 and 12, we showed the effect of Pd on dimensionless moving layer thickness with generalised time Fo. From these figures, we observed that the dimensionless moving layer thickness decreases as the Pd increases for Ki=0.25 at x=1 and 2/3, respectively. That means moving layer thickness increases as the difference of surface and melting temperature and thermal conductivity K of the material increases. Further, the effect of Ki on moving layer thickness λ(Fo) is shown in Figures 13 and 14 for Pd=0.25 at x=1 and x=2/3, respectively. In these figures, we observed that as we increase the value of Ki, the moving layer thickness λ(Fo) decreases as time increases. Similarly, the moving layer thickness increases as the difference of surface and melting temperature and thermal conductivity K of the material increases and decreases as heat flux q increases.

Figure 11:

Effect of Pd on moving layer thickness λ(Fo) with dimensionless time Fo, x=1, Ki=0.25 Case (I).

Figure 12:

Effect of Pd on moving layer thickness λ(Fo) with dimensionless time Fo, x=2/3, Ki=0.25 Case (II).

Figure 13:

Effect of Ki on moving layer thickness λ(Fo) with dimensionless time Fo, x=1, Pd=0.25 Case (I).

Figure 14:

Effect of Ki on moving layer thickness λ(Fo) with dimensionless time Fo, x=2/3, Pd=0.25 Case (II).

In the third kind boundary condition, the effect of Biot number Bi on moving layer thickness λ(Fo) is shown in Figures 15 and 16. In Case (I) at x=1, as we increase the value of Bi, the dimensionless moving layer thickness λ(Fo) increases. The moving layer thickness becomes constant when dimensionless time increases. In Case (II) at x=2/3, as the value of Bi increases, the dimensionless moving layer thickness λ(Fo) increases. By the definition of Bi, as heat transfer coefficient α increases the moving layer thickness increases.

Figure 15:

Effect of Bi on moving layer thickness λ(Fo) with dimensionless time Fo in third kind boundary condition at x=1 Case (I).

Figure 16:

Effect of Bi on moving layer thickness λ(Fo) with dimensionless time Fo in third kind boundary condition at x=2/3 Case (II).

The temperature distribution in first, second and third kind boundary conditions are shown in Figures 17–22, respectively. In these figures, we plotted a graph between temperature distribution with space coordinate x and y for different values of λ obtained at time Fo=0.1. In Case (I) (Figs. 17–19), we observed that as moving layer thickness increases with time the temperatures defined in first, second and third kind boundary conditions increase as space coordinate y increase and the temperature is highest in the third type boundary condition and lowest in the first kind boundary condition. That means the moving layer thickness moves fast in the third kind boundary condition as compared to second and highest in the second kind boundary condition in comparison to the first kind. The material takes less time for complete melting of the third kind boundary condition as compared to the second type boundary condition and less time in the second kind boundary condition in comparison to the first kind boundary condition. In Case (II) (Figs. 20–22), the temperature distribution in first, second and third kind boundary conditions decrease as space coordinate y increases and Figures 20–22 has an opposite behaviour due to the initial conditions.

Figure 17:

Dimensionless temperature θ in first kind boundary condition at λ=0.0041, Fo=0.1 Case (I).

Figure 18:

Dimensionless temperature θ in second kind boundary condition at λ=0.0173, Fo=0.1 Case (I).

Figure 19:

Dimensionless temperature θ in third kind boundary condition at λ=0.0237, Fo=0.1 Case (I).

Figure 20:

Dimensionless temperature θ in first kind boundary condition at λ=0.2018, Fo=0.1 Case (II).

Figure 21:

Dimensionless temperature θ in second kind boundary condition at λ=0.2299, Fo=0.1 Case (II).

Figure 22:

Dimensionless temperature θ in third kind boundary condition at λ=0.3599, Fo=0.1 Case (II).

5 Conclusions

In this present work, the two-dimensional moving boundary problem was solved by Legendre wavelets modified Petrov–Galerkin method. In Petrov–Galerkin discretisation process, we considered the trial function as the linear combination of two-dimensional Legendre wavelets basis function and the test function as two-dimensional Legendre wavelets basis functions for the differential equation. For initial and boundary conditions, we considered the test function as the first element of two-dimensional Legendre wavelets basis function. The convergence analysis of the present method shows that the solution obtained by the present method converges to exact solution for a large value of shifting and dilating parameters. This methodology is extremely efficient to provide the analytical solutions of two-dimensional moving boundary problems. The exceptional accuracy prompts us to conclude that this method is an excellent alternative to other methods for solving two-dimensional moving boundary problems. The effect of parameters such as Pd, Ki and Bi on moving layer thickness λ(Fo) was analysed. Our simulation shows that:

1. The moving layer thickness moves fast in the third kind boundary condition as compared to second and first kind boundary conditions for different values of the particular function g(x, y).

2. The time required for melting of solid is less in boundary condition of the third kind as compared to the boundary condition of the second kind and less in boundary condition of the second kind in comparison to the first kind.

3. As Ki and Pd increase, the moving layer thickness decreases while as Bi increases the moving layer thickness increases.

This method can also be applied to the solution of two-dimensional moving boundary problem which arises during the solidification of the melt.