Nowadays, to investigate the mass, momentum, and heat transfer phenomena, computational fluid dynamics (CFD) study numerical solutions which is associated with the phenomena of thermodynamics, turbo machinery, chemical manufacturing, power generation, weather simulation,  biological enginee-ring, meteorology, aerospace, reaction chemistry, predicting fluid flow and heat transfer (Islam et al., 2019). Though the fundamental source of almost all CFD problems is Navier-Stokes equations which describe many single-phase (liquid or gas, but not both) fluid flows.

Recently, due to the revolution of computer technology, abundant computational grid techniques have been developed which is very efficient to solve numerous engineering problems (Cheniguel and Reghioua, 2013; Chuathong and Toutip, 2011). Among these numerical grid techniques, finite difference method (FDM) is used as a common numerical technique to solve numerous engineering problems (Lau and Kuruganty, 2010). Again, finite element method (FEM) is another kind of commonly used numerical technique which has been applied to solve many heat transfer problems (Patil and Prasad, 2013). A generalized transfer equation for a dependent variable ϕ which can be mass, concentration, heat and momentum is given by Patankar (1980), where discretization technique is applied for CFD analysis (Islam and Hossain, 2019).

Moreover, the finite volume method (FVM) is becoming most popular numerical technique which is commonly used in commercial CFD software such as COMSOL Multiphysics (Versteeg and Malalasekera, 2007). Though to solve any problem these three methods have its own merits besides its demerits, FVM technique is one of the utmost flexible and multipurpose techniques to solve CFD (Uddin et al., 2020). In this paper, we have described an engineering problem by the point of view of FVM. Due to its wide-spread popularity in CFD, this method has been applied previously in different mechanical engineering problem as like thermo elastic and linear elastic problems (Berezovski and Maugin, 2001; Demirdžić and Muzaferija, 1994; and Jasak and Weller, 2000).

The rest of our paper is prepared as follows. A short review of FVM with the help of Tri-Diagonal Matrix Algorithm TDMA is given in Section 2. In Section 3, the numerical solutions obtained by this technique is given where a comparison between exact and our numerical solutions is also described. Finally, in Section 5 we have concluded this paper.

Consider the steady state diffusion of a property T in a one-dimensional domain defined in Fig 1. The process is governed by

d/dx (ψ dT/dx)+q=0                           (1)

Where ψ is the thermal connectivity, T is the temperature and q is the source term (i.e., the rate of heat generation per unit volume). 

2.1 Grid Generation

In the finite volume method first of all, we have to divide the domain into discrete control volumes. Let us consider a number of nodal points in the domain space from A to B. The boundaries (or faces) of control volumes are positioned mid-way between adjacent nodes. Thus each node is surrounded by a control volume or cell. Today it is common practice to set up control volumes near to the edge of the domain in such a way where the physical boundaries coincide with the control volume boundaries. At this point it is appropriate to establish a system of notation that can be used in future developments. 

Fig 1: Grid point generation for finite volume method with control volume.

A general nodal point is identified by O and its neighbors nodes to the north and south are represented by N and S respectively. The south side control volume face is mentioned by s and the north side face of the control volume is mentioned to by n. The distances between the nodes O and S, and between nodes N and O are represented by 〖δx〗_OS and 〖δx〗_NO respectively. In the same way, distances between point O and face s is denoted by 〖δx〗_Os and distances between face n and point O is denoted by〖 δx〗_nO.

2.2 Discretization

To discretize the governing equation, we will integrate the governing equation (1) over the control volume at its nodal point P.

∫_∆V▒〖d/dx (ψ dT/dx)dV+∫_∆V▒〖q dV〗=0〗

Then we get,

(ψA dT/dx)_s-(ψA dT/dx)_n+qA δx=0               (2)

Where A the cross-sectional area of the control volume face is, ∆V is the volume and q is the source term over the control volume. This is a very attractive feature of this finite volume method where the discretized equation contains a clear physical interpretation. Equation 2 states that the diffusive flux of T leaving the south face minus the diffusive flux of T entering the north face is equal to the generation of T over the control volume. In order to derive useful forms of the discretized equations, the interface diffusion coefficient ψ and the gradient dT/(dx )  at west n and east s are required. And the diffusive flux terms are evaluated as -

(ψA dT/dx)_s=ψ_s A((T_S-T_O)/〖δx〗_OS )             (3)

(ψA dT/dx)_n=ψ_n A((T_O-T_N)/〖δx〗_NO )           (4)

For practical situations, the source term q may be a function of the dependent variable. By substituting the values of equations (3) and (4) into equation (2), we get

ψ_s A((T_S-T_O)/〖δx〗_OS )-ψ_n A((T_O-T_N)/〖δx〗_NO )+qA δx=0      (5)

By rearranging equation (5), we can write

((ψ_s  A)/〖δx〗_OS +(ψ_n  A)/〖δx〗_NO ) T_O=((ψ_n  A)/〖δx〗_NO ) T_N+((ψ_s  A)/〖δx〗_OS ) T_S

       +qA δx                                                            (6)

By representing the coefficient of〖 T〗_O, T_N and T_S as〖 a〗_O, a_N and a_(S )respectively then the equation (6) can be written as 

a_O T_O=〖a_N  T〗_N+a_S T_S+S_u                               (7) 

Where (here 〖δx〗_NO=〖δx〗_OS=δx)

a_N a_S a_O S_O S_u

(ψ_n  A)/δx (ψ_s  A)/δx a_N+a_S-S_O 0 qA δx

But there will be a little difference between the boundary conditions at first node (at point A) and for last node (at point B). In this case, linear approximation will be use between the adjacent nodal point and the boundary. For this reasons, we will introduce the linear approximation for temperatures between A and O as like:

ψ_s A((T_S-T_O)/δx)-ψ_A  A((T_O-T_A)/(δx/2))+qA δx=0

⟹((ψ_s  A)/δx+(ψ_A  A)/(δx/2)) T_O=((ψ_s  A)/δx) T_S+((ψ_A  A)/(δx/2)) T_A

                                              +qA δx                     (8) 

This can be rewrite as, 

a_O T_O=〖a_N T〗_N+a_S T_S+S_u                           (9) 

Where 

Similarly, we will introduce the linear approximation for temperatures between O and B as like:

ψ_B  A((T_B-T_O)/(δx/2))-ψ_n A((T_O-T_N)/δx)+qA δx=0   (10) 

⟹((ψ_B  A)/(δx/2)+(ψ_n  A)/δx) T_O=

((ψ_n  A)/δx) T_N+((ψ_B  A)/(δx/2)) T_B+qA δx

This can be rewrite as, 

a_O T_O=〖a_N T〗_N+a_S T_S+S_u                           (11) 

Where 

When we will get the all discretized equations from equation (7), (9) and (11), to solve a problem we have to set up at each of the nodal points. The resulting system of linear algebraic equations will be solved to obtain the distribution of the property T at each nodal point. There are different types of matrix solution technique (Gaussian Elimination Method, Tri-Diagonal Matrix Algorithm, Gauss-Jordan Elimination, etc.) exist which can be used to solve the system of linear algebraic equations. From this system, we will use Gaussian Elimination (GE) Method (Sasaki and Murao, 1982) and Tri-Diagonal Matrix Algorithm (TDMA) (El-Mikkawy, 2004) to find the values of temperature distributions T.

Consider a plane wall whose thickness is 1m (Fig 2) and this has a source effect in this heat conduction system. The two ends of this wall maintained at constant temperature as like 300°C and 800°C respectively. The thermal conductivity ψ of this wall is equals to 0.5 W/m.K, heat generation q = 1000 W/m^3, and the cross-sectional area A is 10 m^2. The main goal of this work is to calculate the temperature distributions of this wall by using our control volume method where we will consider the dimensions of the Y-and Z-directions are so large. For this reason the temperature gradients has significant change along the X -direction only.

Fig 2: A uniform plane wall of 1m length.

3.1 Numerical Solution by FVM

First of all, we have to divide the whole length of the wall into some control volume (along only X-directions). In this work, we will calculate the temperature distribution in five internal points. That is, the difference between two control volumes is 0.2 m (since total length of the wall is 1m, so δx=1/5 m=0.2 m).

Fig 3: Grid point generation of total domain with control volume and half control volume.

For our given problem (one-dimensional heat conduction) governing equation is,

d/dx (ψ dT/dx)+q=0                                            (12)

Where, ψ is the thermal conductivity, T is the temperature and q is the source term. To get the discretization equation, we have to integrate equation (12) over the control volume. After simplifying (which has already described in Section 2) we get as like equation (6),

((ψ_s  A)/〖δx〗_OS +(ψ_n  A)/〖δx〗_NO ) T_O=((ψ_n  A)/〖δx〗_NO ) T_N+((ψ_s  A)/〖δx〗_OS ) T_S

                                       +qA δx                 (13)

By representing the coefficient of T_O, T_N and T_S as a_O, a_N and a_S respectively then the equation (13) can be written as

a_O T_O=〖a_N  T〗_N+a_S T_S+S_u                   (14) 

Where (here 〖δx〗_NO=〖δx〗_OS=δx)

In this problem, the given data are 

By using these values, we get 〖 a〗_N=a_S=25,〖 a〗_O=50, and S_u=2000. Again by using the value of 〖 a〗_N, a_S,〖 a〗_O and S_u  into equation (14), we get the discretized equation for the nodal points 2, 3, and 4 as, 

├ █(50 T_2=25 T_1+25 T_3+2000@50 T_3=25 T_2+25 T_4+2000@50 T_4=25 T_3+25 T_5+2000)}              (15)

On the other hand, at nodal point 1 (at point A) which has already described in Section 2 the coefficients will be,

By using the given values, we get〖 a〗_N=0,〖 a〗_S=25,〖 S〗_O=50, a_O=75, and S_u=17000. Then the equation (14) can be written as, 

75 T_1=25 T_2+17000                              (16) 

Which is the discretized equation for the nodal point 1. Similarly, at nodal point 5 which are a half control volume and the coefficients will be (which has already described in Section 2),

By using the given values, we get of a_N=25,  a_S=0, S_O=50, a_O=75, and S_u=42000. Then the equation (14) can be written as, 

75 T_5=25 T_4+42000                                     (17) 

Which is the discretized equation for the nodal point 5.

Therefore, the set of discretized equations are:

├ █(75 T_1=25 T_2+17000@50 T_2=25 T_1+25 T_3+2000@50 T_3=25 T_2+25 T_4+2000@50 T_4=25 T_3+25 T_5+2000@75 T_5=25 T_4+42000)}           (18)

This system of algebraic linear equation can be solved by different method such as TDMA or Gaussian elimination (GE) method on by hand. But in this we have solved this system of linear equations easily by MATLAB. Then we get our temperature distribution as,

[█(T_1@T_2@T_3@T_4@T_5 )]=[█(450@670@810@860@850)]

3.2 Exact Solution

Our governing equation is  d/dx (ψ dT/dx)+q=0 where boundary conditions, T_A=300°C at x=0 m , and  T_B=800° C,x=1 m. By simplifying this we get our general solution,

T=Cx-(q/2ψ) x^2+D                                   (19)

Where C and D are integral constant. By applying the boundary condition we get C=1500 and D=300, so the complete solution is, 

T=1500x-(q/2ψ) x^2+300                       (20)

For different values of x, we will get the temperature distribution at different nodal points (Table 1).

Fig 4: Graphical representation of (a) Numerical Solution, and (b) Exact Solution.

On the other hand, Fig 4 shows graphical representations of our numerical solution by FVM and exact/analytical solution. Again, Table 2 shows the details presentation of our work where FVM has completed by TDMA and Gaussian elimination (GE) method. Though the two methods show the same result but these numerical solutions have some errors ((Exact Soltution-Numerical Solution)/(Exact Solution)). 

But if we divide the domain into more and more control volume, the error will be minimized than the before.

Table 1: Exact/Analytic solutions for temperature distributions.

Table 2: Comparison table between numerical and exact solutions.

Fig 5: Graphical comparison between control volume method and exact/analytic solution.

In this work, to investigate one dimensional conductive heat transfer in a uniform plane wall, we have used the finite volume method (FVM). Obtaining discretized equation by FVM, we have applied Tri-Diagonal Matrix Algorithm (TDMA) and Gaussian elimination (GE) method. Finally, with the help of MATLAB R2014a, we represent a tabular and a graphical comparison between our numerical solution and exact solutions. The FVM has given outstanding results which have few errors with respect to exact solution but this error can be minimized by taking more and more grid points. That is, this method is very effective, accurate, reliable and easier to appliance in MATLAB or any other programming languages compared to the other exorbitant methods.

We are so much obligated to Md. Sirajul Islam for his patience and thoughtful guidance as a complete supervisor to continue this research work. 

The authors declare that they have no competing interests with respect to the research.