Solution of Unsteady Heat Conduction in an Insulated Rod using the Orthogonal Collocation Method​
​Authors: Housam Binous and Ahmed Bellagi
We compute the first and second derivative matrices using shifted Legendre polynomials (Finlayson, The method of Weighted Residuals and Variational Principles, Academic Press, 1972).
​
The second derivative matrix is then used to solve the unsteady heat conduction in a rod problem (Ramirez, Computational Methods for Process Simulation, 2 nd Edition, Butterworth - Heinemann, 1997).
​
Finally, the temperature profile is compared to the solution of the diffusion equation using NDSolve.
Reference : ​
H. Binous and A. Bellagi, Orthogonal Collocation Methods using MATHEMATICA© in the Graduate Chemical Engineering Curriculum, Computer Applications in Engineering Education Journal, 24 (1), pp. 101 - 113, 2016.
Orthogonal Collocation:​
First and Second Derivative Matrices (A and B) using Shifted Legendre Polynomials (Finlayson, The method of Weighted Residuals and Variational Principles, Academic Press, 1972)
We choose nine internal collocation points (N=9). The code could be changed in order to consider larger numbers of collocation points.
In[]:=
Np=11;
In[]:=
sol=Solve[JacobiP[Np-2,0,0,2x-1.]0,x]
Out[]=
{{x0.0159199},{x0.0819844},{x0.193314},{x0.337873},{x0.5},{x0.662127},{x0.806686},{x0.918016},{x0.98408}}
In[]:=
tbl=Select[Table[x/.sol[[i]],{i,1,Length[sol]}],#>0&];​​x[Np]=1.0;x[1]=10^-20;For[i=2,i<Np,{x[i]=tbl[[i-1]],i++}];
Matrix Q and
-1
Q
In[]:=
Q=Chop@N@Table[Table[x[j]^(i-1),{i,1,Np}],{j,1,Np}]
Out[]=
{{1.,0,0,0,0,0,0,0,0,0,0},{1.,0.0159199,0.000253443,4.03478×
-6
10
,6.42331×
-8
10
,1.02258×
-9
10
,0,0,0,0,0},{1.,0.0819844,0.00672145,0.000551054,0.0000451779,3.70388×
-6
10
,3.03661×
-7
10
,2.48955×
-8
10
,2.04104×
-9
10
,1.67334×
-10
10
,0},{1.,0.193314,0.0373704,0.00722423,0.00139655,0.000269973,0.0000521896,0.000010089,1.95035×
-6
10
,3.7703×
-7
10
,7.28852×
-8
10
},{1.,0.337873,0.114158,0.0385711,0.0130321,0.00440321,0.00148773,0.000502663,0.000169836,0.0000573832,0.0000193882},{1.,0.5,0.25,0.125,0.0625,0.03125,0.015625,0.0078125,0.00390625,0.00195312,0.000976562},{1.,0.662127,0.438412,0.290284,0.192205,0.127264,0.0842649,0.055794,0.0369427,0.0244608,0.0161961},{1.,0.806686,0.650742,0.524944,0.423465,0.341603,0.275566,0.222295,0.179323,0.144657,0.116693},{1.,0.918016,0.842753,0.77366,0.710232,0.652004,0.59855,0.549478,0.504429,0.463074,0.425109},{1.,0.98408,0.968414,0.952997,0.937825,0.922895,0.908203,0.893744,0.879516,0.865514,0.851735},{1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.}}
In[]:=
IQ=Inverse[Q]
Out[]=
{{1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{-90.9989,101.165,-13.3612,4.72783,-2.47227,1.6249,-1.26144,1.13274,-1.19287,1.63597,-0.999604},{2069.97,-2851.38,1052.95,-405.821,217.696,-144.647,112.913,-101.702,107.281,-147.251,89.9893},{-20459.7,30303.3,-14815.9,7688.59,-4474.32,3075.11,-2441.4,2219.44,-2353.19,3237.96,-1979.89},{108568.,-166347.,92670.2,-56972.3,37351.,-27104.3,22129.8,-20432.4,21851.4,-30193.6,18479.4},{-342336.,534442.,-320374.,218661.,-157927.,122257.,-103576.,97695.8,-105754.,147001.,-90088.1},{672661.,-1.06228×
6
10
,666635.,-487656.,379140.,-311918.,275555.,-266815.,293322.,-410893.,252248.},{-832246.,1.32432×
6
10
,-857025.,658140.,-541274.,469506.,-432647.,431482.,-483193.,683350.,-420415.},{630620.,-1.00881×
6
10
,667184.,-530863.,456045.,-413724.,396772.,-408194.,466807.,-667670.,411836.},{-267406.,429426.,-288576.,235847.,-209698.,197565.,-196527.,208587.,-244049.,353618.,-218788.},{48619.2,-78304.1,53262.2,-44443.3,40622.5,-39513.,40622.6,-44443.6,53262.7,-78304.9,48619.7}}
Matrix C
In[]:=
Ci=Chop@N@Table[Table[(i-1)x[j]^(i-2),{i,1,Np}],{j,1,Np}]
Out[]=
{{0,1.,0,0,0,0,0,0,0,0,0},{0,1.,0.0318398,0.000760328,0.0000161391,3.21166×
-7
10
,6.1355×
-9
10
,1.13956×
-10
10
,0,0,0},{0,1.,0.163969,0.0201643,0.00220422,0.000225889,0.0000222233,2.12563×
-6
10
,1.99164×
-7
10
,1.83694×
-8
10
,1.67334×
-9
10
},{0,1.,0.386629,0.112111,0.0288969,0.00698274,0.00161984,0.000365327,0.0000807119,0.0000175531,3.7703×
-6
10
},{0,1.,0.675747,0.342475,0.154284,0.0651607,0.0264193,0.0104141,0.0040213,0.00152853,0.000573832},{0,1.,1.,0.75,0.5,0.3125,0.1875,0.109375,0.0625,0.0351562,0.0195312},{0,1.,1.32425,1.31524,1.16114,0.961024,0.763584,0.589854,0.446352,0.332484,0.244608},{0,1.,1.61337,1.95223,2.09978,2.11732,2.04962,1.92896,1.77836,1.6139,1.44657},{0,1.,1.83603,2.52826,3.09464,3.55116,3.91202,4.18985,4.39582,4.53986,4.63074},{0,1.,1.96816,2.90524,3.81199,4.68913,5.53737,6.35742,7.14995,7.91564,8.65514},{0,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.}}
Matrix A
In[]:=
A=Ci.IQ
Out[]=
{{-90.9989,101.165,-13.3612,4.72783,-2.47227,1.6249,-1.26144,1.13274,-1.19287,1.63597,-0.999604},{-39.0014,30.8986,10.2963,-3.19973,1.61153,-1.04256,0.802941,-0.717885,0.754168,-1.0331,0.631076},{11.1339,-22.2531,5.55383,7.49522,-2.98066,1.77481,-1.31475,1.15149,-1.19623,1.62986,-0.994435},{-5.65889,9.93191,-10.7646,1.96657,6.32292,-2.89898,1.9497,-1.63037,1.65373,-2.22814,1.35615},{3.54228,-5.98715,5.12387,-7.56821,0.724669,5.99958,-3.08403,2.3337,-2.26009,2.983,-1.80762},{-2.4609,4.09377,-3.22466,3.6675,-6.34119,-0.0000185484,6.34123,-3.66755,3.22472,-4.09386,2.46095},{1.80757,-2.98292,2.26003,-2.33366,3.08399,-5.99954,-0.724714,7.56826,-5.12395,5.98727,-3.54235},{-1.3561,2.22805,-1.65367,1.63032,-1.94965,2.89893,-6.32286,-1.96665,10.7647,-9.93208,5.659},{0.994333,-1.62969,1.19611,-1.15139,1.31464,-1.7747,2.98053,-7.49504,-5.55407,22.2534,-11.1342},{-0.630936,1.03287,-0.754008,0.717741,-0.802795,1.0424,-1.61133,3.19948,-10.2959,-30.8992,39.0017},{0.999983,-1.63659,1.19331,-1.13315,1.26186,-1.62539,2.47287,-4.7286,13.3622,-101.167,91.}}
In[]:=
Ci.IQ//MatrixForm
Out[]//MatrixForm=
-90.9989
101.165
-13.3612
4.72783
-2.47227
1.6249
-1.26144
1.13274
-1.19287
1.63597
-0.999604
-39.0014
30.8986
10.2963
-3.19973
1.61153
-1.04256
0.802941
-0.717885
0.754168
-1.0331
0.631076
11.1339
-22.2531
5.55383
7.49522
-2.98066
1.77481
-1.31475
1.15149
-1.19623
1.62986
-0.994435
-5.65889
9.93191
-10.7646
1.96657
6.32292
-2.89898
1.9497
-1.63037
1.65373
-2.22814
1.35615
3.54228
-5.98715
5.12387
-7.56821
0.724669
5.99958
-3.08403
2.3337
-2.26009
2.983
-1.80762
-2.4609
4.09377
-3.22466
3.6675
-6.34119
-0.0000185484
6.34123
-3.66755
3.22472
-4.09386
2.46095
1.80757
-2.98292
2.26003
-2.33366
3.08399
-5.99954
-0.724714
7.56826
-5.12395
5.98727
-3.54235
-1.3561
2.22805
-1.65367
1.63032
-1.94965
2.89893
-6.32286
-1.96665
10.7647
-9.93208
5.659
0.994333
-1.62969
1.19611
-1.15139
1.31464
-1.7747
2.98053
-7.49504
-5.55407
22.2534
-11.1342
-0.630936
1.03287
-0.754008
0.717741
-0.802795
1.0424
-1.61133
3.19948
-10.2959
-30.8992
39.0017
0.999983
-1.63659
1.19331
-1.13315
1.26186
-1.62539
2.47287
-4.7286
13.3622
-101.167
91.
Matrix D
In[]:=
Di=Chop@N@Table[Table[(i-1)(i-2)x[j]^(i-3),{i,1,Np}],{j,1,Np}]
Out[]=
{{0,0,2.,0,0,0,0,0,0,0,0},{0,0,2.,0.0955193,0.00304131,0.0000806955,1.92699×
-6
10
,4.29485×
-8
10
,9.11647×
-10
10
,0,0},{0,0,2.,0.491907,0.0806574,0.0110211,0.00135534,0.000155563,0.000017005,1.79247×
-6
10
,1.83694×
-7
10
},{0,0,2.,1.15989,0.448445,0.144485,0.0418964,0.0113389,0.00292262,0.000726407,0.000175531},{0,0,2.,2.02724,1.3699,0.771421,0.390964,0.184935,0.0833127,0.0361917,0.0152853},{0,0,2.,3.,3.,2.5,1.875,1.3125,0.875,0.5625,0.351562},{0,0,2.,3.97276,5.26094,5.80568,5.76615,5.34509,4.71883,4.01717,3.32484},{0,0,2.,4.84011,7.8089,10.4989,12.7039,14.3473,15.4317,16.0053,16.139},{0,0,2.,5.50809,10.113,15.4732,21.307,27.3842,33.5188,39.5624,45.3986},{0,0,2.,5.90448,11.621,19.0599,28.1348,38.7616,50.8593,64.3496,79.1564},{0,0,2.,6.,12.,20.,30.,42.,56.,72.,90.}}
Matrix B
In[]:=
B=Di.IQ
Out[]=
{{4139.94,-5702.76,2105.9,-811.642,435.392,-289.294,225.825,-203.403,214.561,-294.503,179.979},{2489.47,-3272.99,947.931,-233.77,109.569,-68.7098,52.0827,-46.1577,48.2543,-65.9441,40.2613},{-147.932,426.487,-457.444,217.899,-56.4013,28.2029,-19.1339,15.966,-16.1465,21.714,-13.2103},{36.288,-72.9106,151.043,-206.083,112.347,-30.307,15.9858,-11.7283,11.068,-14.3986,8.69603},{-15.8339,28.5149,-32.6212,93.7385,-140.759,82.7063,-23.4921,13.3381,-11.067,13.5556,-8.07984},{9.84358,-16.9136,15.4284,-23.9171,78.2253,-125.333,78.2253,-23.9171,15.4285,-16.9138,9.84369},{-8.07979,13.5556,-11.067,13.3381,-23.4921,82.7063,-140.759,93.7386,-32.6214,28.5153,-15.8341},{8.69607,-14.3987,11.0681,-11.7284,15.986,-30.3072,112.348,-206.084,151.043,-72.9116,36.2886},{-13.2114,21.716,-16.148,15.9674,-19.1354,28.2046,-56.4036,217.902,-457.448,426.493,-147.937},{40.2731,-65.9634,48.2681,-46.1705,52.0961,-68.7252,109.588,-233.795,947.965,-3273.04,2489.51},{179.997,-294.533,214.583,-203.424,225.847,-289.319,435.424,-811.684,2105.96,-5702.85,4140.}}
In[]:=
Di.IQ//MatrixForm
Out[]//MatrixForm=
4139.94
-5702.76
2105.9
-811.642
435.392
-289.294
225.825
-203.403
214.561
-294.503
179.979
2489.47
-3272.99
947.931
-233.77
109.569
-68.7098
52.0827
-46.1577
48.2543
-65.9441
40.2613
-147.932
426.487
-457.444
217.899
-56.4013
28.2029
-19.1339
15.966
-16.1465
21.714
-13.2103
36.288
-72.9106
151.043
-206.083
112.347
-30.307
15.9858
-11.7283
11.068
-14.3986
8.69603
-15.8339
28.5149
-32.6212
93.7385
-140.759
82.7063
-23.4921
13.3381
-11.067
13.5556
-8.07984
9.84358
-16.9136
15.4284
-23.9171
78.2253
-125.333
78.2253
-23.9171
15.4285
-16.9138
9.84369
-8.07979
13.5556
-11.067
13.3381
-23.4921
82.7063
-140.759
93.7386
-32.6214
28.5153
-15.8341
8.69607
-14.3987
11.0681
-11.7284
15.986
-30.3072
112.348
-206.084
151.043
-72.9116
36.2886
-13.2114
21.716
-16.148
15.9674
-19.1354
28.2046
-56.4036
217.902
-457.448
426.493
-147.937
40.2731
-65.9634
48.2681
-46.1705
52.0961
-68.7252
109.588
-233.795
947.965
-3273.04
2489.51
179.997
-294.533
214.583
-203.424
225.847
-289.319
435.424
-811.684
2105.96
-5702.85
4140.
Unsteady Heat Conduction in an Insulated Rod:​​Must solve following PDE:
∂
t
u=
∂
x,x
u​​ withu(0,t)=u(1,t)=0andu(x,0)=1
​​(Ramirez, Computational Methods for Process Simulation, 2nd Edition, Butterworth-Heinemann,1997).
Temperatures at internal collocation points are
T
i
[t]fori=2to8
In[]:=
DT=Table[D[
T
i
[t],{t,1}],{i,1,Np}]
Out[]=
{
′
T
1
[t],
′
T
2
[t],
′
T
3
[t],
′
T
4
[t],
′
T
5
[t],
′
T
6
[t],
′
T
7
[t],
′
T
8
[t],
′
T
9
[t],
′
T
10
[t],
′
T
11
[t]}
In[]:=
Te=Table[
T
i
[t],{i,1,Np}]
Out[]=
{
T
1
[t],
T
2
[t],
T
3
[t],
T
4
[t],
T
5
[t],
T
6
[t],
T
7
[t],
T
8
[t],
T
9
[t],
T
10
[t],
T
11
[t]}
Writing the system of ODEs to be solved
In[]:=
f[a_,b_]:=ab
In[]:=
eq=MapThread[f,{DT,B.Te}]
In[]:=
eq7=Take[eq,{2,Np-1}]
Initial conditions. Initial temperature is equal to one.
Temperatures at boundary of rod are constant and they are equal to zero.
Equations and initial conditions
Solving the set of ODEs
Temperatures at collocation points at t=0.03
Temperatures profile in the rod at t=0.03
In[]:=
plt1=ListPlot[tbl2,PlotStyle->{Darker@Red,PointSize[0.02]},AspectRatio->1,GridLines->Automatic,Frame->True,ImageSize->400 {1,1},
FrameLabel->{Style["position",14],Style["temperature",14]}]
Solving the PDE using Mathematica
In[]:=
solution =
NDSolve[{D[u[x,t],{t,1}]D[u[x, t],{x,2}],
u[x, 0]1, u[0, t]Exp[-1000 t], u[1, t]Exp[-1000 t]},
u, {x, 0, 1}, {t,0,0.03}]
In[]:=
plt2=Plot[u[x,t]/.solution/.t->0.03,{x,0,1},PlotStyle->{Darker@Blue,PointSize[0.02]},AspectRatio->1,GridLines->Automatic,
Frame->True,ImageSize->400 {1,1},FrameLabel->{Style["position",14],Style["temperature",14]}]
Solutions using Mathematica NDSolve and orthogonal collocation method are identical
In[]:=
Show[plt1,plt2]