Given Information:
Consider, time is represented by variable t.
Length of the rod is 10 cm.
Δx=1 cmΔt=0.1 sT(x=0)=100 °CT(x=10)=50 °C
Here, Δx is the distance between two consecutive nodes along the length of the rod, and Δt is the difference between two consecutive time coordinates.
At, t=0, the temperature of the rod is equivalent to 0°C.
Constant λ defined as k(ΔtΔx2) is equivalent to 0.0835.
Formula used:
At boundary node (0,j), the first derivative in the x dimension is given by finite divided difference as,
∂T∂x≅T1,j−T−1,j2Δx
From Crank Nicolson method, the difference equation for the nodes (except first and the last node) is expressed as,
−λTi−1l+1+2(1+λ)Til+1−λTi+1l+1=λTi−1l+2(1−λ)Til+λTi+1l
Here, i denotes the node varying from 0 to 9, such that i=1 denotes the second node at x=1 cm. Also, l represents the coordinates for time.
Calculation:
For the first step at t=0 s,
Difference equation at node i=1 is solved as,
−λT1−10+1+2(1+λ)T10+1−λT1+10+1=λT1−10+2(1−λ)T10+λT1+10−λT01+2(1+λ)T11−λT21=λT00+2(1−λ)T10+λT20
2(1+λ)T11−λT21=λT00+2(1−λ)T10+λT20+λT01 …… (1)
Substitute λ=0.0835, T10=0, T00=100, T20=0, and T01=100, in equation (1).
2(1+0.0835)T11−(0.0835)T21=(0.0835)(100)+2(1−0.0835)(0)+(0.0835)(0)+(0.0835)(100)
2.167T11−0.0835T21=16.7 …… (2)
Since, T30=0, T40=0, T50=0, T60=0, T70=0, T80=0, T90=0.
Therefore, the difference equations for node i=2,3,4,5,6,7,8, and 9 at t=0 s can be solved as,
−λT11+2(1+λ)T21−λT31=λT10+2(1−λ)T20+λT30
−0.0835T11+2.167T21−0.0835T31=0 …… (3)
−λT21+2(1+λ)T31−λT41=λT20+2(1−λ)T30+λT40
−0.0835T21+2.167T31−0.0835T41=0 …… (4)
−λT31+2(1+λ)T41−λT51=λT30+2(1−λ)T40+λT50
−0.0835T31+2.167T41−0.0835T51=0 …… (5)
−λT41+2(1+λ)T51−λT61=λT40+2(1−λ)T50+λT60
−0.0835T41+2.167T51−0.0835T61=0 …… (6)
−λT51+2(1+λ)T61−λT71=λT50+2(1−λ)T60+λT70
−0.0835T51+2.167T61−0.0835T71=0 …… (7)
−λT61+2(1+λ)T71−λT81=λT60+2(1−λ)T70+λT80
−0.0835T61+2.167T71−0.0835T81=0 …… (8)
−λT71+2(1+λ)T81−λT91=λT70+2(1−λ)T80+λT90
−0.0835T71+2.167T81−0.0835T91=0 …… (9)
−λT81+2(1+λ)T91−λT101=λT80+2(1−λ)T90+λT100−λT81+2(1+λ)T91−λT101=λT80+2(1−λ)T90+λT100−0.0835T81+2.167T91−(0.0835×50)=(0.0835×50)
−0.0835T81+2.167T91=8.35…… (10)
Further, the system of equations (2), (3),(4), (5), (6), (7), (8), and (9) can be written in matrix form as,
[2.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167][T11T21T31T41T51T61T71T81T91]=[16.700000008.35]
Thus,
[T11T21T31T41T51T61T71T81T91]=[2.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167]−1×[16.700000008.35]
Execute the following code in MATLAB to evaluate the results in the above matrix equation
A = [2.167,-0.0835,0,0,0,0,0,0,0;...
-0.0835,2.167,-0.0835,0,0,0,0,0,0;...
0,-0.0835,2.167,-0.0835,0,0,0,0,0;...
0,0,-0.0835,2.167,-0.0835,0,0,0,0;...
0,0,0,-0.0835,2.167,-0.0835,0,0,0;...
0,0,0,0,-0.0835,2.167,-0.0835,0,0;...
0,0,0,0,0,-0.0835,2.167,-0.0835,0;...
0,0,0,0,0,0,-0.0835,2.167,-0.0835;...
0,0,0,0,0,0,0,-0.0835,2.167]
B = [16.7;0;0;0;0;0;0;0;8.35]
X = inv(A)*B
The output values thus obtained are:
X =
7.71798307903248
0.297836314531497
0.011493490904694
0.000443862599904317
2.56738137808539e-05
0.000222426675102992
0.00574676456883212
0.14891815800345
3.85899153954466
Therefore,
[T11T21T31T41T51T61T71T81T91]=[7.7179830.29780.01140.000440.000020.00020.00570.14893.8589]
Furthermore, for the second step to get the values at t=0.2 s, the values of T11,T21,T31,T41,T51,T61,T71,T81, and T91 calculated in first step are multiplied by a factor of 4and incorporated in the original matrix to reflect the values of T calculated in first step. Thus, the system of equations is now written as,
[2.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167][T12T22T32T42T52T62T72T82T92]=4[7.7179830.29780.01140.000440.000020.00020.00570.14893.8589]
Thus, solve as,
[T12T22T32T42T52T62T72T82T92]=[2.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167−0.0835−0.08352.167]−1×[30.87191.19130.045970.0017750.000100.000890.02300.595615.436]
Execute the following code in MATLAB to evaluate the results in the above matrix equation
A = [2.167,-0.0835,0,0,0,0,0,0,0;...
-0.0835,2.167,-0.0835,0,0,0,0,0,0;...
0,-0.0835,2.167,-0.0835,0,0,0,0,0;...
0,0,-0.0835,2.167,-0.0835,0,0,0,0;...
0,0,0,-0.0835,2.167,-0.0835,0,0,0;...
0,0,0,0,-0.0835,2.167,-0.0835,0,0;...
0,0,0,0,0,-0.0835,2.167,-0.0835,0;...
0,0,0,0,0,0,-0.0835,2.167,-0.0835;...
0,0,0,0,0,0,0,-0.0835,2.167]
B =
[30.8719;1.1913;0.04597;0.001775;0.00010;0.00089;0.0230;0.5956; 15.436]
X = inv(A)*B
The output values thus obtained are:
X =
14.2888708601456
1.10279226269952
0.0638337299126952
0.00328788964551248
0.000236412145095749
0.00164989620385863
0.0319231695766006
0.551374157358937
7.14445765673257
Therefore,
[T12T22T32T42T52T62T72T82T92]=[14.28891.10280.06380.00330.000230.001650.03190.55147.1445]