Exercising FVM on a 1D-steady Heat conduction problem using Python

Assume there's a straight rod of 10m length, the left end of which is maintained at a fixed temperature of 100 degree celsius and right end which gives out a constant heat flux of 1W/`m^2` .

The thermal conductivity of the rod varies from left(x=0m) to right(x=10m) as : `k(x)=1+0.01x^2`

Assuming the heat flow a 1-D steady state heat conduction through the rod, Write a computer program to solve problem using Finite Volume Method and Tri-diagonal matrix algorithm

Ans)

Governing eqn. for steady 1-D heat conduction: `frac {del}{del x} (k.frac{del T}{del x})=0`

Exact Method:

`frac {del}{del x} (k.frac{del T}{del x})=0;T_a=100 ; Q_b = 1 ; k(x)=1+0.01x^2`

`=> T(x)=100-10tan^-1(0.1x)`

Finite Volume Method:

Integrating at an inner node P we get,

`int frac {del}{del x} (k.frac{del T}{del x}).dx =0`

`k_e frac{del T_e}{del X} -k_w frac{del T_w}{del X}=0`

By Discretising

`k_e frac{T_E - T _P}{dx} approx k_w frac {T_P-T_W}{dx}`

General form to all nodes hence become

`a_WT_W+a_PT_P+a_ET_E=S`

Note:

at first node : `{del T_w}/{del x} approx {T_P - T_a}/{dx/2}`

at last node : `{del T_e}/{del x} = Q_b`

Hence

A Python Program written to solve the given problem:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

"""
Problem specs:
Governing eqn(Steady 1-D Heat conduction eqn):
d(T)/dx = d/dx(k.dT/dx)

Domain lies along x-axis only

X--->

BC:T=100____________________BC:Q_out=1.0
	(A)_________________________(B)
			K(x)=1+(x/10)^2

Length of domain=10m

Discretised eqn. for eac nodes will be of the form

aw.T_W + ap.T_P + ae.T_E =S
S-source term
aw-west side nodal coefficient
ap-nodal coefficient(self)
ae-east side nodal coefficient

"""
x_start,x_finish=0.0,10.0
lx=x_finish-x_start
nx=10	#no. of finite volumes using to solve the problem
Ta=100.0#left-End(A) fixedTemperature BC
Qb=1.0	#right-End(B) fixedHeatFlux BC
dx=lx/nx#uniform-interval between nodes

#Function for Varying Thermal conductivity 
def K(x):
	return(1+0.01*x**2)

#west-face thermal conductivity of n'th element 
def kw(n):
	return(K(x_start+(n-1)*dx))

#east-face thermal conductivity of n'th element 
def ke(n):
	return(K(x_start+n*dx))

#coefficients of Discretised form
def aw(n):
	if(n==1):
		return(0.0)
	else:
		return(kw(n))
def ap(n):
	if(n==1):
		return(-(2*kw(1)+ke(1)))
	elif(n==nx):
		return(-kw(n))
	else:
		return(-(kw(n)+ke(n)))
def ae(n):
	if(n==nx):
		return(0.0)
	else:
		return(ke(n))
def S(n):
	if(n==1):
		return(-2*kw(1)*Ta)
	elif(n==nx):
		return(Qb*dx)
	else:
		return(0.0)

#forming required Tri-Diagonal matrix 
Tdm=np.array(\
[map(aw,np.arange(1.0,nx+1))\
,map(ap,np.arange(1.0,nx+1))\
,map(ae,np.arange(1.0,nx+1))\
,map(S,np.arange(1.0,nx+1)) ])

T=Tdm[3]*0.0
Tdm=np.transpose(Tdm)
c=Tdm[0][1]
Tdm[0]=Tdm[0]/c
del c

#performing TDMA on the Tri-diagonal matrix
for i in range(1,nx):
	Tdm[i][2]=Tdm[i][2]/(Tdm[i][1]-Tdm[i][0]*Tdm[i-1][2])
	Tdm[i][3]=(Tdm[i][3]-Tdm[i][0]*Tdm[i-1][3])/(Tdm[i][1]-Tdm[i][0]*Tdm[i-1][2])

#performing Back-substitution for solutions
T[nx-1]=Tdm[nx-1][3]

for i in range(nx-2,-1,-1):
	T[i]=Tdm[i][3]-Tdm[i][2]*T[i+1]#FVM solutions

xx=np.linspace(0.0,10.0,nx)#x-axis values for domain length
exacT=100-10*np.arctan(0.1*xx)#exact solution found using Pen&Paper


figure1=plt.plot(xx,T,".",linewidth=0.5,color='red',label='FVM soln.')
figure2=plt.plot(xx,exacT,linewidth=1.5,color='green',label='Exact soln.')
plt.title("Temperature Profile")
plt.xlabel("X")
plt.ylabel("Temp.(in celsius)")

#plt.legend([figure1,figure2],['1','2'])
red_patch = mpatches.Patch(color='red', label='FVM soln.')
#plt.legend(handles=[red_patch])
green_patch = mpatches.Patch(color='green', label='EXACT soln.')
plt.legend(handles=[red_patch,green_patch])

plt.show()

Outputs:

Plot when Cells used in FVM, nx = 10

Plot when Cells used in FVM, nx =50

Note: The last two plots shows the improvement of accuracy with in rease in no. of nodal points