a comparative study of numerical techniques for 2d transient heat ...

Report 8 Downloads 86 Views
International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

A COMPARATIVE STUDY OF NUMERICAL TECHNIQUES FOR 2D TRANSIENT HEAT CONDUCTION EQUATION USING FINITE ELEMENT METHOD By 1

SHARANJEET DHAWAN, 2SHEO KUMAR 1 Department of mathematics 2 Dr. B.R.Ambedkar National Institute of Technology, Jalandhar-144011, India

ABSTRACT Temperature decay in an aluminium plate is observed using Galerkin finite element method for 2D transient heat conduction equation. Taking ∆t of 0.001 hr , temperature variation is studied using unconditionally stable first order and second order accurate schemes - backward Euler and modified CrankNicholson respectively. A comparative study has been made taking different combinations of meshes and numerical schemes. Desired calculations are made taking into account the properties of aluminium plate. Keywords: Temperature, Heat conduction, time, Finite Element Method (FEM). 1.

NOMENCLATURE

T: Temperature t: Time Κ: Thermal conductivity α=

κ

ρC

: Thermal diffusivity

ρ: Density C : Heat capacity Ω: Domain T : Initial temprature ∂Ω : Boundary of the domain 2. INTRODUCTION 2D parabolic equations play a very important role in many applications to study the temperature distribution. Though the analytical solution for the two dimensional heat flow problem is available in the literature but to analyze some phenomena, we need to formulate some model problem correspondingly. There exist many techniques to construct such models and finding their solution. Perhaps no other family of approximation methods has had a greater impact on the theory and practice of numerical methods other than FEM [1]. It has enough solid theoretical foundation that provides more reliability and in many cases makes it possible to analyze mathematically and to estimate errors in approximation. This is why application of Finite Element Method has attracted many linear situations. Because of having different geometries, realistic approximations can be made. Xu and Shen [2] determined transient temperature fields generated by a pulsed laser film in film-substrate considering temperature dependence of the material 38

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

properties. FEM technique and integration for the time variable has been achieved using one order derivative of temperature. Nihad Dukhan [3] presented a one-dimensional heat transfer model for open-cell metal foam, combining the conduction in the ligaments and the convection to the coolant in the pores. Comini and Guide [see 4] proposed a general applicable approach using non-linear physical properties and boundary conditions for transient heat conduction problem using triangular elements for space discretization and using Crank-Nicholson algorithm for each time step. [5] Solved the two dimensional parabolic problem by considering heat conduction in a slab. A space-time finite element has been applied using linear hexahedral elements in space-time domain. Sutradhar [6] found transient temperature distribution for homogeneous and non-homogeneous materials using Laplace transform Galerkin boundary element method. [7] Studied Finite Element Weighted Residual technique for non-linear two-dimensional heat problems using rectangular prism. In the present work a two-dimensional transient heat flow has been considered for solids. A mathematical model has been constructed so that temperature variation can be studied everywhere inside the domain. Solution is started with reformulation of the given differential equation as an equivalent variational problem. The special feature of the finite element method is that the functions are chosen to be piecewise polynomials. Triangular and rectangular finite elements are used. Comparative study has been made taking different combinations of meshes and the appropriate space-time FEM techniques. 3

DESCRIPTION OF NUMERICAL SCHEME

The general partial differential equation for conductive heat flow in two spatial variables describing unsteady temperature distribution in a solid region Ω be defined by 0 x L and 0 y L. Considering a rectangular plate whose edges are bounded by the the line x 0, x L, y 0, y H. Mathematical formulation of the problem is given by the equation ρc

∂T x, y, t ∂t

T x, 0, t

κ

∂ T x, y, t ∂x

T 0, y, t

∂ T x, y, t ∂y

T x, H, t

f x, y, t

T L, y, t , T x, y, 0

T

With the given initial temperature. Using the method of separation of variables, the series solution for this problem is , ,

Where

∑∞

1

∑∞

1

sin

1

sin (1)

1

(2)

To formulate a FEM model of the governing differential equation, we begin with the weak formulation by multiply given equation with the test function v and thus we get ρc

vdΩ

dΩ

κ

fvdΩ

Let 0 x x and 0 y y be the partition of the spatial domain L smaller elements called finite elements for the refinement. We define T x, y, t

N x, y T t

39

(3) H into

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

∂T x, y, t ∂x

∂N x, y T t ∂x

∂T x, y, t ∂y

∂N x, y T t ∂y

Where N x, y linearly independent standard basis are are functions and T t are unknown nodal temperature yet to be determined. Thus we get ∑

T t



dΩ

κ

t

ρcN N dΩ

fN dΩ

(4)

As the weak formulation is a variational statement of the given problem in which we integrated against a test function, the choice of the test function is made in such a way so that after discretization, resulting matrices can be easily solved. We divide the whole domain into a finite number of elements assuming a function for the variation of the variable T in each element. 4.

DISCRETIZATION

In general while discretizing a geometrical domain, the region is divided into quadrilaterals first. The requirement of triangles makes the work quite easy to subdivide the quadrilateral. The solution domain when divided into linear triangles called finite elements, we assume that temperature varies linearly as: T

x, y

C

C

x

C

y

(5)

Thus we approximate unknown nodal temperature using the piecewise approximation: T Where N x, y ; i

x, y, t



T t N x, y

(6)

1,2,3 are the shape functions for the triangle.These shape functions are defined as [8]

Depending upon the situation sometimes we do not feel comfortable with the triangles. Then in the simpler form our quadrilateral becomes a rectangular element with sides parallel to a coordinate system. The interpolation function for a linear rectangular element with sides a and b is as: T x, y C C x C y C xy (7) Following the above procedure adopted for the triangular element, we can proceed for the rectangular element T x, y, t



N x, y T t

(8)

Where N x, y ; i 1,2,3,4 are the shape functions for the rectangular element defined in [8]. Discretizing the whole domain into these non-overlapping finite element and using the shape functions of these elements as nodal points, we get the final assembled form as C

T t

KT t

F

Where {F} is system vector obtained from boundary integration, C K

κ

(9) ρc N NdΩ is mass matrix,

dΩ is stiffness matrix and temperature is assumed to vary linearly over the

time interval ∆t to make approximations. Now to solve the equation for the given initial condition, we use θ family of approximation. where 0 θ 1 is weight factor. When θ 1 it gives backward Euler and for θ 1/2, we get modified Crank-Nicolson.

40

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

Let t t t … t t be the partition of the temporal domain. Defining ∆t t t as the length of the time interval. Starting with the initial condition and solving for the time t , taking θ 1, we have T T θ 1 T . The above discretization process for the time domain and spatial θ θ domain will lead to the system of algebraic equations in the final matrix form as: C

∆t K

u

u

C

∆t F

(10)

Where we can find the unknown nodal temperature at each time level. But this method is only first-order accurate in time. On the other hand sometimes we wish to have more accuracy to take larger time steps, and reducing the round off errors. So, by choosing θ 1/2 we get a second order accurate technique modified Crank Nicolson technique. Assuming the linear variation of the temperature with the time having time interval ∆ and using modified Crank-Nicolson for solving equation (9) we take 2T

T

T

and thus we obtain C

∆t K

T

C

T

∆t F

(11)

Which gives us oscillation free formulation for the linear heat flow equation. Such a formulation is used in various applications which form the basis for further simulations such as laser-induced transient temperature field in film substrate system [2]. 5. RESULTS AND DISCUSSIONS Figure 1 indicates the levels in one symmetric section of the plate for which temperature variation is plotted in forthcoming graphs. Exact solution of the governing differential equation is plotted in Fig.2 showing temperature variation at different levels. It has been observed that there is more cooling (less temperature) at Level 1 and temperature is gradually increasing at center as we move towards Level 5 Similar trends were observed from numerical techniques Figs.3-6. Figures 3-4 represent the results obtained from Backward Euler using rectangular and triangular meshes respectively. Figures 5-6 represent the results obtained from Crank Nicolson (rectangular and triangular mesh). In all the figures maximum temperature was obtained at the center because of cooling from all the sides. Figures 7-8 show the comparison of the temperature variation for different simulation techniques with exact solution. It has also been noticed that numerical schemes are over predicting the exact solution. Figure 9 shows the variation of absolute errors for different techniques plotted at Level 3. Triangular mesh is giving better results as compared to the rectangular mesh for both the differencing techniques. It has also been observed that Backward Euler using triangular mesh is giving minimum error in comparison to others. 6. CONCLUSION Above formulation is quite simple and easy to handle a large variety of linear and non-linear equations. Given domain is discretized into a number of 2D elements. Standard shape functions are used to approximate the nodal temperature at each node inside the domain. Depending upon the requirement, choice of mesh can be selected for the approximation. Proposed idea is flexible enough for different configurations of the material. Testing for the square and rectangular domain gave equally good approximations. Time stepping has been chosen in a way to avoid the divergent oscillations in the temperature. From the results of the case study it is evident that the FEM model formulated is predicting results very close to the exact solution. The work can be extended to some more complex situations for evaluating the temperature distribution in transient heat conduction problems.

41

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

Figure1 Plate configuration showing different solution levels (all the dimensions are in cm).

Level 1 Level 2 Level 3 Level 4 Level 5

25

Temperature(oC)

20

15

10

5

Exact Solution 0 0

50

100

150

200

250

300

Length (cm)

Figure 2 Variation of temperature at different levels (exact solution).

42

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

Level 1 Level 2 Level 3 Level 4 Level 5

30

20

o

Temperature( C)

25

15

10

5

Backward (Rec) 0 0

50

100

150

200

250

300

Length (cm)

Figure 3 Temperature variation at different levels using Backward Euler (Rectangular mesh)

Level 1 Level 2 Level 3 Level 4 Level 5

30

20

o

Temperature( C)

25

15

10

5

Backward(Tri) 0 0

50

100

150

200

250

Length (cm)

Figure 4 Temperature variation at different levels using Backward Euler (Triangular mesh)

43

300

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________ Level 1 Level 2 Level 3 Level 4 Level 5

30

Temperature(oC)

25

20

15

10

5

CrankNicolson(Rec) 0 0

50

100

150

200

250

300

Length (cm)

Figure 5 Temperature variation at different levels using Crank Nicolson (Rectangular mesh)

30

Level 1 Level 2 Level 3 Level 4 Level 5

20

o

Temperature( C)

25

15

10

5

CrankNicolson(Tri) 0 0

50

100

150

200

250

Length (cm)

Figure 6 Temperature variation at different levels using CrankNicolson (Triangular mesh)

44

300

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

25 Exact Solution Crank Nicolson (Rect) CrankNicolson (Tri)

o

Temperature( C)

20

15

10

5

Comparison at Level 3 0 0

50

100

150

200

250

300

Length (cm)

Figure 7 Comparison of temperature with exact solution at Level 3 for CrankNicolson using Rectangular and Triangular mesh

25 Exact Solution Backward (Tri) Backward (Rec)

o

Temperature( C)

20

15

10

5

Comparison at Level 3 0 0

50

100

150

200

250

Length (cm)

Figure 8 Comparison of temperature with exact solution at Level 3 for Backward Euler using Triangular and Rectangular mesh

45

300

International Journal of Research and Reviews in Applied Sciences ISSN: 2076-734X, EISSN: 2076-7366 Volume 1, Issue 1(October 2009) ___________________________________________________________________________________________

3.5

Backward (Tri) CrankNicolson (Tri) Backward (Rec) CrankNicolson (Rec)

3.0

o

Temperature( C)

2.5

2.0

1.5

1.0

Comparison at Level 3

0.5

0.0 0

50

100

150

200

250

300

Length (cm)

Figure 9 Comparison of error obtained by Backward and CrankNicolson using Triangular mesh and Rectangular mesh

6.

REFERENCES

[1] J.Tinsley Oden, “Finite element methods”, Texas Institute for Computatioonal Mechanics, IEEE, 1986. [2] B.Q.Xu, Z.H.Shen, J.Lu, X.W.Ni, S.Y.Zhang, “Numerical simulation of laser - induced transient temperature field in film - substrate system by finite element method”; International Journal of Heat and Mass Transfer, Vol. 46, 2003, pp. 4963-4968. [3] Nihad Dukhan, Pablo D.Quinones - Ramos, Edmundo Cruz - Ruiz, Miguel Velez - Reyes, Elaine P.Scott, “One-dimensional heat transfer analysis in open - cell 10-ppi metal foam”; International Journal of Heat and Mass Transfer, Vol. 48, 2005, pp. 5112-5120. [4] G. Comini and S. Del Guide , “Finite element solution of non-linear heat conduction problems with special reference to phase change”; International Journal for Numerical Methods in Engineering, Vol. 8, 1974, pp. 613-624. [5] Chi-Kyung KIM, “on the numerical computation for solving the two dimensional parabolic equations by space-time finite element method”; JSME International Journal, Vol. 44, 2001, pp. 434-438. [6] Alok Sutradhar, Glaucio H.Paulino, L.J.Gray,”Transient heat conduction in homogeneous and nonhomogeneous materials by Laplace transform Galerkin boundary element method, Engineering”; Analysis with Boundary Elements. Vol. 26, 2002, pp. 119-132. [7] John C.Bruch, JR. and George Zyvoloski, “Transient two dimensional heat conduction problems solved by the finite element method”; International Journal for Numerical Methods in Engineering, Vol. 8, 1974 pp. 481-494. [8] O.C.Zienkiewicz, R.L.Taylor & J.Z.Zhu, “The Finite element Method Its Basis & Fundamentals”. Elsevier, Sixth edition.

46