Chapter 10 - Two-Dimensional Steady & Transient Heat Conduction

Exa10.1 Page 195

In [1]:
from __future__ import division
tnew=[]
e=[]
for i in range(1,10):
    tnew.append(101)
    e.append(1)                       #assumed values

t=1e-6
#since all the nodes are interior nodes so discretized eqn used is eqn 10.10
while e[0]>t and e[1]>t and e[2]>t and e[3]>t and e[4]>t and e[5]>t and e[6]>t and e[7]>t and e[8]>t:
    told=[]
    for i in range(1,10):
        told.append(tnew[i-1])
    tnew[0]=(told[1]+40+told[3])/4                       #on solving eqns for various nodes we get,
    tnew[1]=(tnew[0]+told[2]+told[4]+20)/4
    tnew[2]=(tnew[1]+told[5]+420)/4
    tnew[3]=(told[4]+tnew[0]+told[6]+20)/4
    tnew[4]=(tnew[1]+told[7]+told[5]+tnew[3])/4
    tnew[5]=(tnew[2]+tnew[4]+told[8]+400)/4
    tnew[6]=(tnew[3]+told[7]+40)/4
    tnew[7]=(tnew[4]+tnew[6]+told[8]+20)/4
    tnew[8]=(tnew[5]+420+tnew[7])/4
    e=[]
    for i in range(1,10):
        e.append(abs(tnew[i-1]-told[i-1]))
    

print "the values of T from 1st element to last is:\n"
for i in range(1,10):
    print tnew[i-1],'\t',
the values of T from 1st element to last is:

47.142854591 	91.2499974482 	182.857141581 	57.3214260196 	114.999997448 	220.178570153 	47.1428558669 	91.2499987241 	182.857142219 	

Exa10.2 Page 198

In [2]:
from __future__ import division
tnew=[]
e=[]
for i in range(1,6):
    tnew.append(101)
    e.append(1)

t=1e-6
#since there is no source term so we get the following eqns.
while e[0]>t and e[1]>t and e[2]>t and e[3]>t and e[4]>t:
    told=[]
    for i in range(1,6):
        told.append(tnew[i-1])
        
    tnew[0]=(told[1]*2+300)/4
    tnew[1]=(tnew[0]+told[2]+300)/4
    tnew[2]=(tnew[1]+told[3]+200)/4
    tnew[3]=(told[4]+tnew[2]+300)/4
    tnew[4]=(2*tnew[3]+300)/4
    
    for i in range(1,6):
        e[i-1]=(abs(tnew[0]-told[i-1]))
    

print "the values of T from 1st element to last is:\n"
for i in range(1,6):
    print tnew[i-1],'\t',
the values of T from 1st element to last is:

145.833332213 	141.666666106 	120.833333053 	141.666666527 	145.833333263 	

Exa10.3 Page 201

In [3]:
from __future__ import division
from numpy import zeros
tn=zeros([5,5])
for j in range(1,6):
    tn[j-1,0]=400                       #defining conditions

for j in range(2,5):
    tn[0,j-1]=20
    tn[4,j-1]=20
    tn[j-1,4]=20

e=[]    
for i in range(1,10):
    e.append(1)

for i in range(2,5):
    tn[i-1,1:4]=60


t1=1e-6
t=zeros([4,4])
while e[0]>t1 and e[1]>t1 and e[2]>t1 and e[3]>t1 and e[4]>t1 and e[5]>t1 and e[6]>t1 and e[7]>t1 and e[8]>t1:
    for i in range(1,4):
        for j in range(1,4):
            t[i,j]=tn[i,j]
        
    #using red-black gauss-seidel method
    for i in range(1,4):
        for j in range(1,4):
            tn[i,j]=(tn[i+1,j]+tn[i-1,j]+tn[i,j+1]+tn[i,j-1])/4
    #now getting the absolute difference of the 9 variables
    e[0]=abs(t[1,1]-tn[1,1])
    e[1]=abs(t[1,2]-tn[1,2])
    e[2]=abs(t[1,3]-tn[1,3])
    e[3]=abs(t[2,1]-tn[2,1])
    e[4]=abs(t[2,2]-tn[2,2])
    e[5]=abs(t[2,3]-tn[2,3])
    e[6]=abs(t[3,1]-tn[3,1])
    e[7]=abs(t[3,2]-tn[3,2])
    e[8]=abs(t[3,3]-tn[3,3])

#now defining positions of the various nodes according to red black combination
R1=t[1,3]
R2=t[3,3]
R3=t[2,2]
R4=t[1,2]
R5=t[1,3]
B1=t[3,2]
B2=t[2,1]
B3=t[2,3]
B4=t[1,2]
print "the value of RED points respectively : %.4f, %.4f, %.4f, %.4f & %.4f"%(R1,R2,R3,R4,R5)

print "the value of BLUE points respectively : %.4f, %.4f, %.4f & %.4f"%(B1, B2, B3, B4)
the value of RED points respectively : 47.1429, 47.1429, 115.0000, 91.2500 & 47.1429
the value of BLUE points respectively : 91.2500, 220.1786, 57.3214 & 91.2500

Exa10.8 Page 202

In [4]:
from __future__ import division
from numpy import zeros
tnew=[0];told=[0];e=[0];
for i in range(1,10):
    tnew.append(101)
    told.append(0)
    e.append(1)                    #assumed values

t=1e-6
while e[(1)]>t and e[(2)]>t and e[(3)]>t and e[(4)]>t and e[(5)]>t and e[(6)]>t and e[(7)]>t and e[(8)]>t and e[(9)]>t:
    for i in range(1,10):
        told[i]=tnew[(i)]
    #using eqn 10.10 for the interior nodes and convective boundary conditions for corner nodes
    tnew[(1)]=(told[(2)]+50+.5*told[(4)]+100/3)*3/7
    tnew[(2)]=(tnew[(1)]+told[(3)]+told[(5)]+100)/4
    tnew[(3)]=(tnew[(2)]+told[(6)]+600)/4
    tnew[(4)]=(told[(5)]+.5*tnew[(1)]+.5*told[(7)]+100/3)*3/7
    tnew[(5)]=(tnew[(2)]+told[(8)]+told[(6)]+tnew[(4)])/4
    tnew[(6)]=(tnew[(3)]+tnew[(5)]+told[(9)]+500)/4
    tnew[(7)]=(.5*tnew[(4)]+.5*told[(8)]+100/3)*3/4
    tnew[(8)]=(tnew[(5)]+.5*tnew[(7)]+.5*told[(9)]+100/3)*3/7
    tnew[(9)]=(tnew[(6)]+100/3+.5*tnew[(8)]+250)*3/7
    for i in range(1,10):
        e[i]=abs(tnew[i]-told[i])
    

print "the values of T from 1st element to last is"
for i in range(1,10):
    print tnew[i]
the values of T from 1st element to last is
157.831321254
192.469876146
280.72288988
184.939753978
231.325296989
330.421684639
175.903610664
217.469876356
309.638552636

Exa10.9 Page 211

In [5]:
from __future__ import division
from numpy import zeros
tnew=[0];told=[0];e=[0];
for i in range(1,10):
    tnew.append(101)
    told.append(0)
    e.append(1)                    #assumed values

t=1e-6
while e[(1)]>t and e[(2)]>t and e[(3)]>t and e[(4)]>t and e[(5)]>t and e[(6)]>t and e[(7)]>t and e[(8)]>t and e[(9)]>t:
    for i in range(1,10):
        told[(i)]=tnew[(i)]
    #using the basic discretization eqn. for all the nodes since the boundary conditions vary for each point
    tnew[(1)]=(told[(2)]+1.25+told[(4)])/2.05
    tnew[(2)]=(.5*tnew[(1)]+.5*told[(3)]+told[(5)]+1.25)/2.05
    tnew[(3)]=(.5*tnew[(2)]+.5*told[(6)]+1.25)/1.05
    tnew[(4)]=(told[(5)]+.5*tnew[(1)]+45)/2
    tnew[(5)]=(tnew[(2)]+told[(8)]+90+tnew[(4)])/4
    tnew[(6)]=(.5*tnew[(3)]+tnew[(5)]+.5*told[(7)]+91.25)/3.05
    tnew[(7)]=(.5*tnew[(6)]+.5*told[(8)]+91.25)/2.05
    tnew[(8)]=(91.25+.5*tnew[(7)]+.5*told[(9)])/2.05
    tnew[(9)]=(47.125+.5*tnew[(8)])/1.025
    for i in range(1,10):
        e[i]=abs(tnew[i]-told[i])
    

print "the values of T from 1st element to last is"
for i in range(1,10):
    print told[i]
the values of T from 1st element to last is
83.6250164077
83.3426002934
81.8505046198
86.8385228347
86.864496399
86.0434332286
86.7854379748
87.2768602399
88.5496879219

Exa10.10 Page 213

In [6]:
from __future__ import division
from numpy import zeros
print "the soln of eg 10.10-->Alternating Direction Implicit Method"
nmax=25
a=[0,1,1]                                #defining variables
b=[-4,-4,-4]
c=[1,1]
t=zeros([6,6])
tstar=zeros([6,6])
t[1,2]=20 ; t[1,3]=20 ; t[1,4]=20; t[2,1]=20; t[3,1]=20 ; t[4,1]=20 ; t[5,2]=20 ; t[5,3]=20 ; t[5,4]=20 ; t[2,5]=400; t[3,5]=400; t[4,5]=400
tstar[1,2]=20 ; tstar[1,3]=20; tstar[1,4]=20 ; tstar[2,1]=20 ; tstar[3,1]=20 ; tstar[4,1]=20 ; tstar[5,2]=20 ; tstar[5,3]=20 ; tstar[5,4]=20 ; tstar[2,5]=400 ; tstar[3,5]=400 ; tstar[4,5]=400
for i in range(2,5):
    for j in range(2,5):
        t[i,j]=20
    

#solving by TDMA Method
d=zeros(4)
for nn in range(1,nmax+1):
    for jj in range(2,5):
        d[(1)]=-t[1,jj]-t[2,jj+1]-tstar[2,jj-1]
        d[(2)]=-t[3,jj+1]-tstar[3,jj-1]
        d[(3)]=-t[5,jj]-t[4,jj+1]-tstar[4,jj-1]
        i=1 ; n=3
    Beta=[0,b[i]]
    Gamma=[0,d[i]/Beta[i]]
    i1=i+1
    for j in range(i1,n+1):
        Beta.append(b[(j-1)]-a[(j-1)]*c[(j-2)]/Beta[(j-1)])
        Gamma.append((d[(j)]-a[(j-1)]*Gamma[(j-1)])/Beta[(j)])
    #end 
    x=zeros(n)
    x[-1]=Gamma[n-1]
    n1=n-i
    for k in range(1,n1+1):
        j=n-k
        x[j-1]=Gamma[(j)]-c[(j-1)]*x[j-1]/Beta[j]
    #end
    tstar[2,jj]=x[0]
    tstar[3,jj]=x[1]
    tstar[4,jj]=x[2]
#end
    for ii in range(2,5):
        d[(1)]=-t[ii,1]-tstar[ii+1,2]-tstar[ii-1,2]
        d[(2)]=-tstar[ii+1,3]-tstar[ii-1,3]
        d[(3)]=-t[ii,5]-tstar[ii+1,4]-tstar[ii-1,4]
    i=1 ; n=3
    Beta[(i)]=b[i]
    Gamma[i]=d[(i)]/Beta[(i)]
    i1=i+1
    for j in range(i1,n+1):
        Beta[(j)]=b[(j-1)]-a[(j-1)]*c[j-2]/Beta[(j-1)]
        Gamma[(j)]=(d[(j)]-a[(j-1)]*Gamma[(j-1)])/Beta[(j)]
    #end 
    
    x[-1]=Gamma[(n)]
    n1=n-i
    for k in range(1,n1+1):
        j=n-k
        x[j]=Gamma[(j)]-c[(j-1)]*x[(j)]/Beta[(j)]
    #end
    t[ii,2]=x[0]
    t[ii,3]=x[1]
    t[ii,4]=x[2]


print "the soln by ADI method is"
print t[2,4],t[2,3],t[2,2]
print "--------------"
print t[3,4],t[3,3],t[3,2]
print "--------------"
print t[4,4],t[4,3],t[4,2]
the soln of eg 10.10-->Alternating Direction Implicit Method
the soln by ADI method is
20.0 20.0 20.0
--------------
20.0 20.0 20.0
--------------
48.1904761905 43.6666666667 105.0

Exa 10.11 Page 215

In [7]:
from __future__ import division
from numpy import arange,zeros
a=[0,0];b=[0,0];c=[0,0]
for k in range(2,11):
    a.append(-1)
    b.append(2.4)
    c.append(-1)

alpha=1 ; delta_t=.05 ; delta_x=.1
m=alpha*delta_t/delta_x**2
b[1]=2.4
c[1]=-2
t=zeros([12,12])
tstar=zeros([12,12])

for k in range(1,12):
    t[11,k]=400 ; tstar[11,k]=400 ;  t[k,11]=400 ; tstar[k,11]=400


for i in range(1,11):
    for j in range(1,11):
        t[i,j]=0
        tstar[i,j]=0
    
d=zeros(11)
for tm in arange(.05,5+.05,.05):
    for jj in range(1,11):
        if jj==1:
            for ii in range(1,11):
                d[(ii)]=2*t[ii,2]-1.6*t[ii,1]
            d[(10)]=d[(10)]+400
        else:
            for ii in range(1,11):
                d[(ii)]=t[ii,jj+1]+t[ii,jj-1]-1.6*t[ii,jj]
            d[(10)]=d[(10)]+400
        i=1 ; n=10
        Beta=[0,b[i]]
        Gamma=[0,d[i]/Beta[i]]
        i1=i+1
        for j in range(i1,n+1):
            Beta.append(b[(j)]-a[(j)]*c[(j-1)]/Beta[(j-1)])
            Gamma.append((d[(j)]-a[(j)]*Gamma[(j-1)])/Beta[(j)])
    
        x=zeros(n+1)
        x[(n)]=Gamma[(n)] 
        x[(j)]=Gamma[(j)]-c[(j)]*x[(j)]/Beta[(j)]
        for count in range(1,11):
            tstar[count,jj]=x[(count)]
    
    for ii in range(1,11):
        if ii==1:
            for jj in range(1,11):
                d[(jj)]=2*tstar[ii+1,1]-1.6*tstar[ii,1]
            d[(10)]=d[(10)]+400
        else:
            for jj in range(1,11):
                d[(jj)]=tstar[ii-1,jj]+tstar[ii+1,jj]-1.6*tstar[ii,jj]
                d[(10)]=d[(10)]+400
                
        i=1 ;  n=10
        Beta[(i)]=b[i]
        Gamma[(i)]=d[(i)]/Beta[(i)]
        i1=i+1
        for j in range(i1,n+1):
            Beta[(j)]=b[(j)]-a[(j)]*c[(j-1)]/Beta[(j-1)]
            Gamma[(j)]=(d[(j)]-a[(j)]*Gamma[(j-1)])/Beta[(j)]
        
        x[(n)]=Gamma[(n)]
        n1=n-i
        for k in range(1,n1+1):
            j =n-k
            x[(j)]=Gamma[(j)]-c[(j)]*x[(j+1)]/Beta[(j)]
    
        for count in range(1,11):
            t[ii,count]=x[count]
  

print "the soln by ADI is"
for i in range(1,11):
    print t[i,i]
the soln by ADI is
1.58563151796
1.90275782155
2.98098725376
5.25161158747
9.62288055617
17.8433017473
33.2010436374
61.8392029825
945.175268307
-467.149084346