- Email: [email protected]
A finite-difference method of high-order accuracy for the solution of transient nonlinear diffusive–convective problem in three dimensions
Case Studies in Thermal Engineering 3 (2014) 43–50
Contents lists available at ScienceDirect
Case Studies in Thermal Engineering journal homepage: www.elsevier.com/locate/csite
A finite-difference method of high-order accuracy for the solution of transient nonlinear diffusive–convective problem in three dimensions Marco Donisete de Campos a,b,n, Estaner Claro Romão c, Luiz Felipe Mendes de Moura b a Federal University of Mato Grosso, Institute of Exact and Earth Sciences, Av. Governador Jaime Campos, 6390, 78600-000 Barra do Garças, Mato Grosso, Brazil b Mechanical Engineering Faculty, State University of Campinas, PO Box 6122, 13083-970 Campinas, São Paulo, Brazil c Department of Basic and Environmental Sciences, EEL-USP, University of São Paulo, 12602-810 Lorena, São Paulo, Brazil
a r t i c l e i n f o
abstract
Article history: Received 18 November 2013 Received in revised form 13 February 2014 Accepted 7 March 2014 Available online 3 April 2014
This paper presents an efficient technique of linearization of the nonlinear convective terms present in engineering problems involving heat and mass transfer and fluid mechanics. From two numerical applications, this technique with a method of highorder finite differences is validated by numerical solution of transient nonlinear diffusive– convective problem in three dimensions. & 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).
Keywords: Linearization technique Nonlinear problems Finite differences
1. Introduction Numerical simulation of engineering problems has been successfully performed since decades ago. Problems governed by linear differential equations have been solved easily and successfully by most numerical methods proposed in the literature, among them we can mention the works involving, e.g., finite difference method, approximate methods for fractional ordinary differential equations, finite element method, spectral method and matrix approach. Based on the high-order finite difference method, that is the focus of the work, several authors have proposed various methodologies for numerically solving the transient diffusive–convective problems. Among these papers [1], established an exponential high-order compact alternating direction implicit method for the numerical solution of unsteady 2D convection–diffusion problems using the Crank–Nicolson scheme for the time discretization and an exponential fourthorder compact difference formula for the spatial discretization. The unconditionally stable character of the scheme has been verified by a discrete Fourier analysis and the computational results showed the accuracy, efficiency and robustness of the method [2] developed a high order compact alternating direction implicit method for solving two-dimensional unsteady convection–diffusion problems. The method is unconditionally stable and second-order accurate in time and fourth-order accurate in space. Three numerical studies are carried out to demonstrate its high accuracy and efficiency and, although it
n Corresponding author at: Federal University of Mato Grosso, Institute of Exact and Earth Sciences, Av. Governador Jaime Campos, 6390, 78600-000 Barra do Garças, Mato Grosso, Brazil. Tel.: þ 55 663 405 317. E-mail address: [email protected] (M.D. de Campos).
http://dx.doi.org/10.1016/j.csite.2014.03.001 2214-157X/& 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).
44
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
has been proposed for convection–diffusion problems, can also be used to solve pure diffusion or pure convection problems [3] using the compact difference scheme with the operator-splitting technique solved the multidimensional time fractional diffusion problem, investigating the stability and accuracy of the scheme [4] developed an exponential high order compact alternating direction implicit method for solving three dimensional unsteady convection–diffusion equations. The method is fourth-order in space, second order in time, and proved to be unconditionally stable. Three numerical experiments are performed to demonstrate its high accuracy and efficiency and to show its superiority over the classical Douglas–Gunn alternating direction implicit scheme and the Karaa's high-order alternating direction implicit scheme [5]. In this work, the aim is to present a formulation that uses an accurate linearization technique with easy to implement and that uses only one iteration at each time step, optimizing the calculation of the velocity profile. 2. Formulation In this work, we propose a solution, by the high-order finite difference method for the transient nonlinear convection– diffusion equation in three dimensions which is given by 2 ∂u ∂u ∂u ∂u ∂ u ∂2 u ∂2 u þu þu þ u ¼ υ þ 2þ 2 ð1Þ 2 ∂t ∂x ∂y ∂z ∂x ∂y ∂z where u(x,y,z,t) is the velocity field in the x,y,z-directions and υ is the kinematic viscosity. For the temporal discretization, considering the equation given by Eq. (1), we use the method called α family of approximation [6], in which a weighted average of the time derivative of a dependent variable is approximated on two consecutive time steps by linear interpolation of the values of the variable at two steps n n þ 1 ~ n þ 1 fug ~ n ∂u fug ∂u ∂u ¼ ð1 αÞ þ α ð2Þ ∂t ∂t ∂t Δt n þ 1 where 0 rα r1; tA [tn, tn þ 1], n ¼0,1,2,…,mt, where mt is the number of steps in time, and { }n refers to the value of the enclosed quantity at time n and Δt n þ 1 ¼ t n þ 1 t n is the (nþ 1)th time step. Adopting α ¼ 0:5; we obtain the well-known Crank–Nicolson method to carry out the time discretization as follows: nþ1 2 nþ1 u un ∂ u ∂ 2 un þ 1 ∂2 un þ 1 ∂un þ 1 ∂un þ 1 ∂un þ 1 ¼ 0:5 υ un þ 1 un þ 1 þυ þυ un þ 1 2 2 2 Δt ∂x ∂y ∂z ∂x ∂y ∂z 2 n n n n ∂ u ∂2 un ∂2 un n ∂u n ∂u n ∂u u u ð3Þ þ 0:5 υ þυ þυ u ∂x ∂y ∂z ∂x2 ∂y2 ∂z2 In order to linearize the convective terms of the previous equation will be used the technique proposed by [7] according to which for a sufficiently small time step, these terms can be expanded into a Taylor series after the first-derivative terms. The result is as follows: un þ 1
∂un þ 1 ∂un þ 1 ∂un ∂un un þ un þ 1 un : ∂x ∂x ∂x ∂x
ð4aÞ
un þ 1
∂un þ 1 ∂un þ 1 ∂un ∂un un þ un þ 1 un ∂y ∂y ∂y ∂y
ð4bÞ
un þ 1
∂un þ 1 ∂un þ 1 ∂un ∂un un þ un þ 1 un ∂z ∂z ∂z ∂z
ð4cÞ
This method is commonly known as Newton's method since it provides a quadratic convergence [8]. Note that this technique does not require an iterative linearization in each time step, making quicker the calculation of u. And, substituting Eqs. (4a)–(4c) into the convective terms of Eq. (2), it yields nþ1 2 nþ1 u un ∂ u ∂ 2 un þ 1 ∂2 un þ 1 ∂un þ 1 ∂un ∂un ∂un þ 1 ¼ 0:5 υ un þ 1 þun un þυ þυ un Δt ∂x ∂x ∂x ∂y ∂x2 ∂y2 ∂z2 n n nþ1 n n ∂u ∂u ∂u ∂u ∂u un þ 1 þ un un un þ 1 þun ∂y ∂y ∂z ∂z ∂z 2 n n n ∂ u ∂2 un ∂2 un ∂u ∂u ∂un þ 0:5 υ un un þυ þυ un 2 2 2 ∂x ∂y ∂z ∂x ∂y ∂z n n nþ1 n un þ 1 n þ 1 ∂u n þ 1 ∂u n ∂u n þ 1 ∂u þu þu þu þu þ ¼F ð5Þ ∂x ∂y ∂z ∂z Δt where F¼
2 n un ∂ u ∂2 un ∂2 un þ 0:5 υ : þυ þυ 2 2 2 Δt ∂x ∂y ∂z
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
45
Now, discretizing in the space by Central Difference Method of O(Δx2) for nodes with dx, dy or dz distance from the boundary, we have ! ! 1 þ1 1 þ1 uni þþ1;j;k 2uni;j;k þ uni þ1;j;k uni;jþþ11;k 2uni;j;k þuni;jþ11;k 0:5υ 0:5υ Δx2 Δy2 ! ! ! þ1 þ1 þ1 1 1 uni;j;k 2uni;j;k þ uni;j;k uni þþ1;j;k uni þ1;j;k uni;jþþ11;k uni;jþ11;k þ1 1 n n þ 0:5u þ0:5u 0:5υ i;j;k i;j;k 2Δx 2Δy Δz2 ! ! nþ1 nþ1 n n n u u ∂u ∂u ∂u 1 i;j;k þ 1 i;j;k 1 i;j;k i;j;k i;j;k þ1 þ 0:5 þ þ0:5 þ 0:5 uni;j;k ¼ F1 þ 0:5uni;j;k Δt 2Δz ∂x ∂y ∂z ! ! ! n n n 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 ui;j;k 1 þ 2 ui;j 1 þ 2 ui 1;j;k ) 2 Δz Δy Δx Δz Δy Δx ! ! ∂uni;j;k ∂uni;j;k ∂uni;j;k υ υ υ 1 þ1 þ 0:5 þ þ þ þ þ þ uni;j;k ∂x ∂y ∂z Δx2 Δy2 Δz2 Δt ! ! ! n n n 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 ui þ 1;j;k þ 2 þ ui;j þ 1;k þ 2 þ ui;j;k þ 1 ¼ F 1 : þ 2þ ð6Þ Δx Δy Δz Δx Δy Δz where F1 ¼
uni;j;k Δt
þ 0:5υ
uniþ 1;j;k 2uni;j;k þuni 1;j;k Δx2
þ
uni;j þ 1;k 2uni;j;k þ uni;j 1;k Δy2
þ
uni;j;k þ 1 2uni;j;k þ uni;j;k 1 Δz2
! :
Now, considering the internal nodes using the Central Difference Method with O(Δx4), we have 1 1 þ1 1 1 uni þþ2;j;k þ 16uni þþ1;j;k 30uni;j;k þ16uni þ1;j;k uni þ2;j;k
0:5υ þ þ þ
12Δx2
þ1 þ1 þ1 þ1 þ1 uni;j;k þ16uni;j;k 30uni;j;k þ 16uni;j;k uni;j;k þ2 þ1 1 2
12Δy 1 þ0; 5 Δt
∂uni;j;k ∂x
þ
∂uni;j;k ∂y
þ
∂uni;j;k ∂z
þ
þ1 uni;jþþ12;k þ16uni;jþþ11;k 30uni;j;k þ 16uni;jþ11;k uni;jþ12;k
12Δy2
!
12Δz2 uni;jþþ12;k þ8uni;jþþ11;k 8uni;jþ11;k þ uni;jþ12;k
þ
þ0:5uni;j;k
1 1 1 1 uni þþ2;j;k þ 8uni þþ1;j;k 8uni þ1;j;k þuni þ2;j;k
þ1 þ1 þ1 þ1 uni;j;k þ8uni;j;k 8uni;j;k þ uni;j;k þ2 þ1 1 2
12Δx !
12Δz
!! þ1 ¼ F2 uni;j;k
! ! ! uni;j;k uni;j;k n þ 1 uni;j;k υ 2υ υ nþ1 u u un þ 1 ) þ þ þ þ 24Δz2 24Δz i;j;k 2 3Δz2 3Δz i;j;k 1 24Δy2 24Δy i;j 2;k ! ! ! uni;j;k uni;j;k uni;j;k n þ 1 2υ υ 2υ nþ1 nþ1 þ þ þ þ u u u 3Δy2 3Δy i;j 1;k 24Δx2 24Δx i 2;j;k 3Δx2 3Δx i 1;j;k ! ! ! uni;j;k uni;j;k n þ 1 uni;j;k υ 2υ υ nþ1 þ) þ þ þ þ u u un þ 1 24Δz2 24Δz i;j;k 2 3Δz2 3Δz i;j;k 1 24Δy2 24Δy i;j 2;k ! ! ! uni;j;k uni;j;k uni;j;k n þ 1 2υ υ 2υ nþ1 nþ1 þ þ þ þ u u u 3Δy2 3Δy i;j 1;k 24Δx2 24Δx i 2;j;k 3Δx2 3Δx i 1;j;k ! ! ∂uni;j;k ∂uni;j;k ∂uni;j;k 1:25υ 1:25υ 1:25υ 1 þ1 þ 0:5 þ þ þ þ þ þ uni;j;k ∂x ∂y ∂z Δx2 Δy2 Δz2 Δt ! ! ! uni;j;k n þ 1 uni;j;k uni;j;k 2υ υ 2υ nþ1 u u un þ 1 þ þ þ þ þ 3Δx2 3Δx i þ 1;j;k 24Δx2 24Δx i þ 2;j;k 3Δy2 3Δy i;j þ 1;k ! ! ! uni;j;k uni;j;k n þ 1 uni;j;k υ 2υ υ nþ1 þ þ þ þ u u un þ 1 ¼ F 2 24Δy2 24Δy i;j þ 2;k 3Δz2 3Δz i;j;k þ 1 24Δz2 24Δz i;j;k þ 2
ð7Þ
where F2 ¼
uni;j;k Δt þ
þ 0:5υ
uniþ 2;j;k þ16uniþ 1;j;k 30uni;j;k þ 16uni 1;j;k uni 2;j;k 12Δx2
uni;j;k 2 þ 16uni;j;k 1 30uni;j;k þ 16uni;j;k þ 1 uni;j;k þ 2 12Δz2
! :
þ
uni;j 2;k þ 16uni;j 1;k 30uni;j;k þ 16uni;j þ 1;k uni;j þ 2;k 12Δy2
46
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
3. Numerical applications Here, we present two numerical applications: first, from the proposal to an exact solution to Eq. (1), will be considered the greatest error (L1 error) for different spatial and temporal refinements, besides the computational time and the convergence achieved for meshes adopted. All computations were run on a Intel Core i7/2.4G private computer using double precision arithmetic. In the second application, using the sinusoidal velocity profiles, it was demonstrated situations in which softens the mesh refinement of the solution numerical oscillations of the velocity field.
3.1. Application 1 Considering the exact solution of Eq. (1) given by uðx; y; z; tÞ ¼ expððx þ y 2z þ6tÞ=νÞ, was done the comparison with the numerical solution, with ν¼1, t ¼0.1, 0 rx,y,zr0.1 and using a regular mesh with h¼Δx ¼Δy¼Δz¼ 0.1. Table 1 shows the errors associated with the numerical solution compared with the analytical solution according to the L1 norm, considering the maximum error for stopping criterion for the Gauss–Seidel on the order of 10 10, the CPU time and the rate of convergence (CR) defined by CR ¼ ðlogðE1 =E2 Þ= logðΔt 1 =Δt 2 ÞÞ, where E1 and E2 are the L1 norm errors with the temporal grid sizes Δt 1 and Δt 2 , respectively. The boundary and initial conditions are directly taken from this analytical solution. At first we set h¼0.005, t¼0.1 and change Δt from 0.005 to 0.001. The L1 norm errors and the CPU time used for the present method are depicted in Table 1. We notice that the present method achieves the accuracy of 2.03E 08 with a temporal refinement of Δt¼0.001 whereas for h¼ 0.0025 achieves the accuracy of 9.03E 09. Fig. 1 shows the standard error L1 and velocity profile u in the xy-plane, considering Δx ¼Δy¼Δz¼Lx/20 and Lx ¼ Ly ¼Lz ¼0.1. It can be noted that the greatest error occurs in the center of the plane analyzed.
Table 1 L1 error, CPU time and the rate of with t ¼0.1. Δt
0.00500 0.00250 0.00167 0.00125 0.00100
h ¼ 0.005
h ¼ 0.0025
L1 norm
CPU time
CR
L1 norm
CPU time
CR
4.69E 07 1.16E 07 5.29E 08 3.06E 08 2.03E 08
60.781 103.172 137.375 165.328 189.110
– 2.02 1.93 1.91 1.84
4.58E 07 1.05E 07 4.14E 08 1.92E 08 9.03E 09
1744.875 2873.437 3708.547 4842.578 5391.891
– 2.12 2.29 2.67 3.38
Fig. 1. L1 error (left) and profile u (right) both in the xy-plane, for z¼ 0.05 and t¼ 0.1.
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
47
Fig. 2. Two dimensional velocity profile in the xy-plane for z ¼0.5 (left), and in the xz-plane for y¼0.5 (right) considering υ ¼1.
Fig. 3. Two dimensional velocity profile in the xy-plane for z ¼0.5 (left), and in the xz-plane for y¼0.5 (right), considering υ ¼0.1.
3.2. Application 2 Now, considering the computational domain as being an unit cube and t ¼0.1, with the initial condition uðx; y; z; 0Þ ¼ sin ð2πxÞ cos ð2πyÞz and the boundary condition uðx; y; z; tÞ ¼ sin ð2πxÞ cos ð2πyÞzet . Fig. 2 shows the bidimensional profile, respectively, of u-velocity in the xy-plane with z¼0.5 and in the xz-plane with y¼ 0.5 considering υ¼ 1. Figs. 3 and 4 show the two dimensional velocity profile u in the xy-plane with z¼0.5 and in the xz-plane with y¼0.5 are given, respectively, υ ¼0.1 and υ¼0.01, with Δx¼Δy¼Δz ¼Lx/20 and Δt¼ Lt/20. Aiming to decrease the numerical oscillations in the previous cases, the mesh was refined, considering now Δx ¼Δy¼Δz¼Lx/40 and Δt ¼Lt/20. Thus, in the conditions of the previous case, was obtained the velocity profile shown in Fig. 5, in which we note that the refinement softened the numerical oscillations.
48
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
Fig. 4. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right) considering υ ¼ 0.01.
Fig. 5. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right), considering υ ¼0.01.
Figs. 6 and 7 illustrates the numerical results for two dimensional velocity profile u in the xy-plane with z ¼0.5 and in the xz-plane with y¼0.5 are given, respectively, υ ¼0.001 with Δx¼ Δy ¼Δz¼Lx/190 and Δt¼Lt/100 and υ ¼0.0001, with Δx ¼Δy¼Δz¼ Lx/20 and Δt¼Lt/1000. 4. Conclusions In this work, an efficient technique of linearization of the nonlinear convective terms present in engineering problems involving heat and mass transfer and fluid mechanics is presented. This technique is very simple, stable and accurate and technique does not require any type of iterative process within each time step, which represents a major saving on the computational time. It used acopled with the method with a method
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
49
Fig. 6. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right), considering υ ¼ 0.001.
Fig. 7. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right), considering υ ¼0.0001.
of high-order finite differences is validated by two numerical applications of tridimensional transient nonlinear diffusionconvection equation in three dimensions.
Acknowledgements We acknowledged the support from the National Council of Scientific Development and Technology, CNPq, Brazil (Proc. 500382/2011-5) and Mato Grosso Research Foundation, Fapemat, Brazil (Proc. 292470/2010). References [1] Tian ZF, Ge YB. A fourth-order compact ADI method for solving two-dimensional unsteady convection–diffusion problems. J Comput Appl Math 2007;198:268–86. [2] Tian ZF. A rational high-order compact ADI method for unsteady convection–diffusion equations. Comput Phys Commun 2011;182:649–62.
50
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
[3] Cui M. Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation. Numer Algorithms 2013;62:383–409. [4] Ge YB, Tian ZF, Zhang J. An exponential high order compact ADI (EHOC ADI) method for 3D unsteady convection–diffusion problems. Numer Methods Partial Differ Equ 2013;29:186–205. [5] Karaa S. A high-order compact ADI method for solving three-dimensional unsteady convection–diffusion problems. Numer Methods Partial Differ Equ 2006;22:983–93. [6] Reddy JN. An introduction to the finite element method. 3rd ed. McGraw-Hill; 2013. [7] Dennis Jr JE, Schmabel RB. Numerical methods for unconstrained optimization and nonlinear equations. Prentice-Hall; 1983. [8] Jiang BN. The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics. Springer; 1998.
Contents lists available at ScienceDirect
Case Studies in Thermal Engineering journal homepage: www.elsevier.com/locate/csite
A finite-difference method of high-order accuracy for the solution of transient nonlinear diffusive–convective problem in three dimensions Marco Donisete de Campos a,b,n, Estaner Claro Romão c, Luiz Felipe Mendes de Moura b a Federal University of Mato Grosso, Institute of Exact and Earth Sciences, Av. Governador Jaime Campos, 6390, 78600-000 Barra do Garças, Mato Grosso, Brazil b Mechanical Engineering Faculty, State University of Campinas, PO Box 6122, 13083-970 Campinas, São Paulo, Brazil c Department of Basic and Environmental Sciences, EEL-USP, University of São Paulo, 12602-810 Lorena, São Paulo, Brazil
a r t i c l e i n f o
abstract
Article history: Received 18 November 2013 Received in revised form 13 February 2014 Accepted 7 March 2014 Available online 3 April 2014
This paper presents an efficient technique of linearization of the nonlinear convective terms present in engineering problems involving heat and mass transfer and fluid mechanics. From two numerical applications, this technique with a method of highorder finite differences is validated by numerical solution of transient nonlinear diffusive– convective problem in three dimensions. & 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).
Keywords: Linearization technique Nonlinear problems Finite differences
1. Introduction Numerical simulation of engineering problems has been successfully performed since decades ago. Problems governed by linear differential equations have been solved easily and successfully by most numerical methods proposed in the literature, among them we can mention the works involving, e.g., finite difference method, approximate methods for fractional ordinary differential equations, finite element method, spectral method and matrix approach. Based on the high-order finite difference method, that is the focus of the work, several authors have proposed various methodologies for numerically solving the transient diffusive–convective problems. Among these papers [1], established an exponential high-order compact alternating direction implicit method for the numerical solution of unsteady 2D convection–diffusion problems using the Crank–Nicolson scheme for the time discretization and an exponential fourthorder compact difference formula for the spatial discretization. The unconditionally stable character of the scheme has been verified by a discrete Fourier analysis and the computational results showed the accuracy, efficiency and robustness of the method [2] developed a high order compact alternating direction implicit method for solving two-dimensional unsteady convection–diffusion problems. The method is unconditionally stable and second-order accurate in time and fourth-order accurate in space. Three numerical studies are carried out to demonstrate its high accuracy and efficiency and, although it
n Corresponding author at: Federal University of Mato Grosso, Institute of Exact and Earth Sciences, Av. Governador Jaime Campos, 6390, 78600-000 Barra do Garças, Mato Grosso, Brazil. Tel.: þ 55 663 405 317. E-mail address: [email protected] (M.D. de Campos).
http://dx.doi.org/10.1016/j.csite.2014.03.001 2214-157X/& 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).
44
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
has been proposed for convection–diffusion problems, can also be used to solve pure diffusion or pure convection problems [3] using the compact difference scheme with the operator-splitting technique solved the multidimensional time fractional diffusion problem, investigating the stability and accuracy of the scheme [4] developed an exponential high order compact alternating direction implicit method for solving three dimensional unsteady convection–diffusion equations. The method is fourth-order in space, second order in time, and proved to be unconditionally stable. Three numerical experiments are performed to demonstrate its high accuracy and efficiency and to show its superiority over the classical Douglas–Gunn alternating direction implicit scheme and the Karaa's high-order alternating direction implicit scheme [5]. In this work, the aim is to present a formulation that uses an accurate linearization technique with easy to implement and that uses only one iteration at each time step, optimizing the calculation of the velocity profile. 2. Formulation In this work, we propose a solution, by the high-order finite difference method for the transient nonlinear convection– diffusion equation in three dimensions which is given by 2 ∂u ∂u ∂u ∂u ∂ u ∂2 u ∂2 u þu þu þ u ¼ υ þ 2þ 2 ð1Þ 2 ∂t ∂x ∂y ∂z ∂x ∂y ∂z where u(x,y,z,t) is the velocity field in the x,y,z-directions and υ is the kinematic viscosity. For the temporal discretization, considering the equation given by Eq. (1), we use the method called α family of approximation [6], in which a weighted average of the time derivative of a dependent variable is approximated on two consecutive time steps by linear interpolation of the values of the variable at two steps n n þ 1 ~ n þ 1 fug ~ n ∂u fug ∂u ∂u ¼ ð1 αÞ þ α ð2Þ ∂t ∂t ∂t Δt n þ 1 where 0 rα r1; tA [tn, tn þ 1], n ¼0,1,2,…,mt, where mt is the number of steps in time, and { }n refers to the value of the enclosed quantity at time n and Δt n þ 1 ¼ t n þ 1 t n is the (nþ 1)th time step. Adopting α ¼ 0:5; we obtain the well-known Crank–Nicolson method to carry out the time discretization as follows: nþ1 2 nþ1 u un ∂ u ∂ 2 un þ 1 ∂2 un þ 1 ∂un þ 1 ∂un þ 1 ∂un þ 1 ¼ 0:5 υ un þ 1 un þ 1 þυ þυ un þ 1 2 2 2 Δt ∂x ∂y ∂z ∂x ∂y ∂z 2 n n n n ∂ u ∂2 un ∂2 un n ∂u n ∂u n ∂u u u ð3Þ þ 0:5 υ þυ þυ u ∂x ∂y ∂z ∂x2 ∂y2 ∂z2 In order to linearize the convective terms of the previous equation will be used the technique proposed by [7] according to which for a sufficiently small time step, these terms can be expanded into a Taylor series after the first-derivative terms. The result is as follows: un þ 1
∂un þ 1 ∂un þ 1 ∂un ∂un un þ un þ 1 un : ∂x ∂x ∂x ∂x
ð4aÞ
un þ 1
∂un þ 1 ∂un þ 1 ∂un ∂un un þ un þ 1 un ∂y ∂y ∂y ∂y
ð4bÞ
un þ 1
∂un þ 1 ∂un þ 1 ∂un ∂un un þ un þ 1 un ∂z ∂z ∂z ∂z
ð4cÞ
This method is commonly known as Newton's method since it provides a quadratic convergence [8]. Note that this technique does not require an iterative linearization in each time step, making quicker the calculation of u. And, substituting Eqs. (4a)–(4c) into the convective terms of Eq. (2), it yields nþ1 2 nþ1 u un ∂ u ∂ 2 un þ 1 ∂2 un þ 1 ∂un þ 1 ∂un ∂un ∂un þ 1 ¼ 0:5 υ un þ 1 þun un þυ þυ un Δt ∂x ∂x ∂x ∂y ∂x2 ∂y2 ∂z2 n n nþ1 n n ∂u ∂u ∂u ∂u ∂u un þ 1 þ un un un þ 1 þun ∂y ∂y ∂z ∂z ∂z 2 n n n ∂ u ∂2 un ∂2 un ∂u ∂u ∂un þ 0:5 υ un un þυ þυ un 2 2 2 ∂x ∂y ∂z ∂x ∂y ∂z n n nþ1 n un þ 1 n þ 1 ∂u n þ 1 ∂u n ∂u n þ 1 ∂u þu þu þu þu þ ¼F ð5Þ ∂x ∂y ∂z ∂z Δt where F¼
2 n un ∂ u ∂2 un ∂2 un þ 0:5 υ : þυ þυ 2 2 2 Δt ∂x ∂y ∂z
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
45
Now, discretizing in the space by Central Difference Method of O(Δx2) for nodes with dx, dy or dz distance from the boundary, we have ! ! 1 þ1 1 þ1 uni þþ1;j;k 2uni;j;k þ uni þ1;j;k uni;jþþ11;k 2uni;j;k þuni;jþ11;k 0:5υ 0:5υ Δx2 Δy2 ! ! ! þ1 þ1 þ1 1 1 uni;j;k 2uni;j;k þ uni;j;k uni þþ1;j;k uni þ1;j;k uni;jþþ11;k uni;jþ11;k þ1 1 n n þ 0:5u þ0:5u 0:5υ i;j;k i;j;k 2Δx 2Δy Δz2 ! ! nþ1 nþ1 n n n u u ∂u ∂u ∂u 1 i;j;k þ 1 i;j;k 1 i;j;k i;j;k i;j;k þ1 þ 0:5 þ þ0:5 þ 0:5 uni;j;k ¼ F1 þ 0:5uni;j;k Δt 2Δz ∂x ∂y ∂z ! ! ! n n n 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 ui;j;k 1 þ 2 ui;j 1 þ 2 ui 1;j;k ) 2 Δz Δy Δx Δz Δy Δx ! ! ∂uni;j;k ∂uni;j;k ∂uni;j;k υ υ υ 1 þ1 þ 0:5 þ þ þ þ þ þ uni;j;k ∂x ∂y ∂z Δx2 Δy2 Δz2 Δt ! ! ! n n n 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 0:5υ 0:25ui;j;k n þ 1 ui þ 1;j;k þ 2 þ ui;j þ 1;k þ 2 þ ui;j;k þ 1 ¼ F 1 : þ 2þ ð6Þ Δx Δy Δz Δx Δy Δz where F1 ¼
uni;j;k Δt
þ 0:5υ
uniþ 1;j;k 2uni;j;k þuni 1;j;k Δx2
þ
uni;j þ 1;k 2uni;j;k þ uni;j 1;k Δy2
þ
uni;j;k þ 1 2uni;j;k þ uni;j;k 1 Δz2
! :
Now, considering the internal nodes using the Central Difference Method with O(Δx4), we have 1 1 þ1 1 1 uni þþ2;j;k þ 16uni þþ1;j;k 30uni;j;k þ16uni þ1;j;k uni þ2;j;k
0:5υ þ þ þ
12Δx2
þ1 þ1 þ1 þ1 þ1 uni;j;k þ16uni;j;k 30uni;j;k þ 16uni;j;k uni;j;k þ2 þ1 1 2
12Δy 1 þ0; 5 Δt
∂uni;j;k ∂x
þ
∂uni;j;k ∂y
þ
∂uni;j;k ∂z
þ
þ1 uni;jþþ12;k þ16uni;jþþ11;k 30uni;j;k þ 16uni;jþ11;k uni;jþ12;k
12Δy2
!
12Δz2 uni;jþþ12;k þ8uni;jþþ11;k 8uni;jþ11;k þ uni;jþ12;k
þ
þ0:5uni;j;k
1 1 1 1 uni þþ2;j;k þ 8uni þþ1;j;k 8uni þ1;j;k þuni þ2;j;k
þ1 þ1 þ1 þ1 uni;j;k þ8uni;j;k 8uni;j;k þ uni;j;k þ2 þ1 1 2
12Δx !
12Δz
!! þ1 ¼ F2 uni;j;k
! ! ! uni;j;k uni;j;k n þ 1 uni;j;k υ 2υ υ nþ1 u u un þ 1 ) þ þ þ þ 24Δz2 24Δz i;j;k 2 3Δz2 3Δz i;j;k 1 24Δy2 24Δy i;j 2;k ! ! ! uni;j;k uni;j;k uni;j;k n þ 1 2υ υ 2υ nþ1 nþ1 þ þ þ þ u u u 3Δy2 3Δy i;j 1;k 24Δx2 24Δx i 2;j;k 3Δx2 3Δx i 1;j;k ! ! ! uni;j;k uni;j;k n þ 1 uni;j;k υ 2υ υ nþ1 þ) þ þ þ þ u u un þ 1 24Δz2 24Δz i;j;k 2 3Δz2 3Δz i;j;k 1 24Δy2 24Δy i;j 2;k ! ! ! uni;j;k uni;j;k uni;j;k n þ 1 2υ υ 2υ nþ1 nþ1 þ þ þ þ u u u 3Δy2 3Δy i;j 1;k 24Δx2 24Δx i 2;j;k 3Δx2 3Δx i 1;j;k ! ! ∂uni;j;k ∂uni;j;k ∂uni;j;k 1:25υ 1:25υ 1:25υ 1 þ1 þ 0:5 þ þ þ þ þ þ uni;j;k ∂x ∂y ∂z Δx2 Δy2 Δz2 Δt ! ! ! uni;j;k n þ 1 uni;j;k uni;j;k 2υ υ 2υ nþ1 u u un þ 1 þ þ þ þ þ 3Δx2 3Δx i þ 1;j;k 24Δx2 24Δx i þ 2;j;k 3Δy2 3Δy i;j þ 1;k ! ! ! uni;j;k uni;j;k n þ 1 uni;j;k υ 2υ υ nþ1 þ þ þ þ u u un þ 1 ¼ F 2 24Δy2 24Δy i;j þ 2;k 3Δz2 3Δz i;j;k þ 1 24Δz2 24Δz i;j;k þ 2
ð7Þ
where F2 ¼
uni;j;k Δt þ
þ 0:5υ
uniþ 2;j;k þ16uniþ 1;j;k 30uni;j;k þ 16uni 1;j;k uni 2;j;k 12Δx2
uni;j;k 2 þ 16uni;j;k 1 30uni;j;k þ 16uni;j;k þ 1 uni;j;k þ 2 12Δz2
! :
þ
uni;j 2;k þ 16uni;j 1;k 30uni;j;k þ 16uni;j þ 1;k uni;j þ 2;k 12Δy2
46
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
3. Numerical applications Here, we present two numerical applications: first, from the proposal to an exact solution to Eq. (1), will be considered the greatest error (L1 error) for different spatial and temporal refinements, besides the computational time and the convergence achieved for meshes adopted. All computations were run on a Intel Core i7/2.4G private computer using double precision arithmetic. In the second application, using the sinusoidal velocity profiles, it was demonstrated situations in which softens the mesh refinement of the solution numerical oscillations of the velocity field.
3.1. Application 1 Considering the exact solution of Eq. (1) given by uðx; y; z; tÞ ¼ expððx þ y 2z þ6tÞ=νÞ, was done the comparison with the numerical solution, with ν¼1, t ¼0.1, 0 rx,y,zr0.1 and using a regular mesh with h¼Δx ¼Δy¼Δz¼ 0.1. Table 1 shows the errors associated with the numerical solution compared with the analytical solution according to the L1 norm, considering the maximum error for stopping criterion for the Gauss–Seidel on the order of 10 10, the CPU time and the rate of convergence (CR) defined by CR ¼ ðlogðE1 =E2 Þ= logðΔt 1 =Δt 2 ÞÞ, where E1 and E2 are the L1 norm errors with the temporal grid sizes Δt 1 and Δt 2 , respectively. The boundary and initial conditions are directly taken from this analytical solution. At first we set h¼0.005, t¼0.1 and change Δt from 0.005 to 0.001. The L1 norm errors and the CPU time used for the present method are depicted in Table 1. We notice that the present method achieves the accuracy of 2.03E 08 with a temporal refinement of Δt¼0.001 whereas for h¼ 0.0025 achieves the accuracy of 9.03E 09. Fig. 1 shows the standard error L1 and velocity profile u in the xy-plane, considering Δx ¼Δy¼Δz¼Lx/20 and Lx ¼ Ly ¼Lz ¼0.1. It can be noted that the greatest error occurs in the center of the plane analyzed.
Table 1 L1 error, CPU time and the rate of with t ¼0.1. Δt
0.00500 0.00250 0.00167 0.00125 0.00100
h ¼ 0.005
h ¼ 0.0025
L1 norm
CPU time
CR
L1 norm
CPU time
CR
4.69E 07 1.16E 07 5.29E 08 3.06E 08 2.03E 08
60.781 103.172 137.375 165.328 189.110
– 2.02 1.93 1.91 1.84
4.58E 07 1.05E 07 4.14E 08 1.92E 08 9.03E 09
1744.875 2873.437 3708.547 4842.578 5391.891
– 2.12 2.29 2.67 3.38
Fig. 1. L1 error (left) and profile u (right) both in the xy-plane, for z¼ 0.05 and t¼ 0.1.
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
47
Fig. 2. Two dimensional velocity profile in the xy-plane for z ¼0.5 (left), and in the xz-plane for y¼0.5 (right) considering υ ¼1.
Fig. 3. Two dimensional velocity profile in the xy-plane for z ¼0.5 (left), and in the xz-plane for y¼0.5 (right), considering υ ¼0.1.
3.2. Application 2 Now, considering the computational domain as being an unit cube and t ¼0.1, with the initial condition uðx; y; z; 0Þ ¼ sin ð2πxÞ cos ð2πyÞz and the boundary condition uðx; y; z; tÞ ¼ sin ð2πxÞ cos ð2πyÞzet . Fig. 2 shows the bidimensional profile, respectively, of u-velocity in the xy-plane with z¼0.5 and in the xz-plane with y¼ 0.5 considering υ¼ 1. Figs. 3 and 4 show the two dimensional velocity profile u in the xy-plane with z¼0.5 and in the xz-plane with y¼0.5 are given, respectively, υ ¼0.1 and υ¼0.01, with Δx¼Δy¼Δz ¼Lx/20 and Δt¼ Lt/20. Aiming to decrease the numerical oscillations in the previous cases, the mesh was refined, considering now Δx ¼Δy¼Δz¼Lx/40 and Δt ¼Lt/20. Thus, in the conditions of the previous case, was obtained the velocity profile shown in Fig. 5, in which we note that the refinement softened the numerical oscillations.
48
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
Fig. 4. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right) considering υ ¼ 0.01.
Fig. 5. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right), considering υ ¼0.01.
Figs. 6 and 7 illustrates the numerical results for two dimensional velocity profile u in the xy-plane with z ¼0.5 and in the xz-plane with y¼0.5 are given, respectively, υ ¼0.001 with Δx¼ Δy ¼Δz¼Lx/190 and Δt¼Lt/100 and υ ¼0.0001, with Δx ¼Δy¼Δz¼ Lx/20 and Δt¼Lt/1000. 4. Conclusions In this work, an efficient technique of linearization of the nonlinear convective terms present in engineering problems involving heat and mass transfer and fluid mechanics is presented. This technique is very simple, stable and accurate and technique does not require any type of iterative process within each time step, which represents a major saving on the computational time. It used acopled with the method with a method
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
49
Fig. 6. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right), considering υ ¼ 0.001.
Fig. 7. Two dimensional velocity profile in the xy-plane for z¼ 0.5 (left), and in the xz-plane for y¼ 0.5 (right), considering υ ¼0.0001.
of high-order finite differences is validated by two numerical applications of tridimensional transient nonlinear diffusionconvection equation in three dimensions.
Acknowledgements We acknowledged the support from the National Council of Scientific Development and Technology, CNPq, Brazil (Proc. 500382/2011-5) and Mato Grosso Research Foundation, Fapemat, Brazil (Proc. 292470/2010). References [1] Tian ZF, Ge YB. A fourth-order compact ADI method for solving two-dimensional unsteady convection–diffusion problems. J Comput Appl Math 2007;198:268–86. [2] Tian ZF. A rational high-order compact ADI method for unsteady convection–diffusion equations. Comput Phys Commun 2011;182:649–62.
50
M.D. de Campos et al. / Case Studies in Thermal Engineering 3 (2014) 43–50
[3] Cui M. Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation. Numer Algorithms 2013;62:383–409. [4] Ge YB, Tian ZF, Zhang J. An exponential high order compact ADI (EHOC ADI) method for 3D unsteady convection–diffusion problems. Numer Methods Partial Differ Equ 2013;29:186–205. [5] Karaa S. A high-order compact ADI method for solving three-dimensional unsteady convection–diffusion problems. Numer Methods Partial Differ Equ 2006;22:983–93. [6] Reddy JN. An introduction to the finite element method. 3rd ed. McGraw-Hill; 2013. [7] Dennis Jr JE, Schmabel RB. Numerical methods for unconstrained optimization and nonlinear equations. Prentice-Hall; 1983. [8] Jiang BN. The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics. Springer; 1998.