Chapter 30, Problem 5P

Repeat Example 30.3, but for $\Delta x=1\text{ cm}$.

To determine

To calculate: The temperature distribution of a long, thin aluminum rod with the given properties (refer Example 30.3, in the text book) when the rod is initially at $0°\text{C}$ and the derivative boundary conditions at $x=10\text{ cm}$ is 0.

$$

Answer to Problem 5P

Solution: The temperature distribution is,

$\left[\begin{array}{c}{T}_1^2\\ T_2^2\\ T_3^2\\ T_4^2\\ T_5^2\\ T_6^2\\ T_7^2\\ T_8^2\\ T_9^2\end{array}\right]=\left[\begin{array}{c}14.2889\\ 1.1028\\ 0.0638\\ 0.0033\\ 0.00023\\ 0.00165\\ 0.0319\\ 0.5514\\ 7.1445\end{array}\right]$

Explanation of Solution

Given Information:

Consider, time is represented by variable t.

Length of the rod is 10 cm.

$\begin{array}{l}\Delta x=1\text{ cm}\\ \Delta t=0.1\text{ s}\\ \text{T}\left(x=0\right)=100°\text{C}\\ \text{T}\left(x=10\right)=50°\text{C}\end{array}$

Here, $\Delta x$ is the distance between two consecutive nodes along the length of the rod, and $\Delta t$ is the difference between two consecutive time coordinates.

At, $t=0$, the temperature of the rod is equivalent to $0°\text{C}$.

Constant $\lambda$ defined as $k\left(\frac{\Delta t}{\Delta x^2}\right)$ is equivalent to $0.0835$.

Formula used:

At boundary node $\left(0,j\right)$, the first derivative in the x dimension is given by finite divided difference as,

$\frac{\partial T}{\partial x}\cong \frac{T_{1,j}-T_{-1,j}}{2\Delta x}$

From Crank Nicolson method, the difference equation for the nodes (except first and the last node) is expressed as,

$-\lambda T_{i-1}^{l+1}+2\left(1+\lambda \right)T_i^{l+1}-\lambda T_{i+1}^{l+1}=\lambda T_{i-1}^l+2\left(1-\lambda \right)T_i^l+\lambda T_{i+1}^l$

Here, $i$ denotes the node varying from 0 to 9, such that $i=1$ denotes the second node at $x=1\text{ cm}$. Also, $l$ represents the coordinates for time.

Calculation:

For the first step at $t=0\text{ s}$,

Difference equation at node $i=1$ is solved as,

$\begin{array}{c}-\lambda T_{1-1}^{0+1}+2\left(1+\lambda \right)T_1^{0+1}-\lambda T_{1+1}^{0+1}=\lambda T_{1-1}^0+2\left(1-\lambda \right)T_1^0+\lambda T_{1+1}^0\\ -\lambda T_0^1+2\left(1+\lambda \right)T_1^1-\lambda T_2^1=\lambda T_0^0+2\left(1-\lambda \right)T_1^0+\lambda T_2^0\end{array}$

$2\left(1+\lambda \right)T_1^1-\lambda T_2^1=\lambda T_0^0+2\left(1-\lambda \right)T_1^0+\lambda T_2^0+\lambda T_0^1$ …… (1)

Substitute $\lambda =0.0835$, $T_1^0=0$, $T_0^0=100$, $T_2^0=0$, and $T_0^1=100$, in equation (1).

$2\left(1+0.0835\right)T_1^1-\left(0.0835\right)T_2^1=\left(0.0835\right)\left(100\right)+2\left(1-0.0835\right)\left(0\right)+\left(0.0835\right)\left(0\right)+\left(0.0835\right)\left(100\right)$

$2.167T_1^1-0.0835T_2^1=16.7$ …… (2)

Since, $T_3^0=0,T_4^0=0,T_5^0=0,T_6^0=0,T_7^0=0,T_8^0=0,T_9^0=0$.

Therefore, the difference equations for node $i=2,3,4,5,6,7,8,\text{ and 9}$ at $t=0\text{ s}$ can be solved as,

$-\lambda T_1^1+2\left(1+\lambda \right)T_2^1-\lambda T_3^1=\lambda T_1^0+2\left(1-\lambda \right)T_2^0+\lambda T_3^0$

$-0.0835T_1^1+2.167T_2^1-0.0835T_3^1=0$ …… (3)

$-\lambda T_2^1+2\left(1+\lambda \right)T_3^1-\lambda T_4^1=\lambda T_2^0+2\left(1-\lambda \right)T_3^0+\lambda T_4^0$

$-0.0835T_2^1+2.167T_3^1-0.0835T_4^1=0$ …… (4)

$-\lambda T_3^1+2\left(1+\lambda \right)T_4^1-\lambda T_5^1=\lambda T_3^0+2\left(1-\lambda \right)T_4^0+\lambda T_5^0$

$-0.0835T_3^1+2.167T_4^1-0.0835T_5^1=0$ …… (5)

$-\lambda T_4^1+2\left(1+\lambda \right)T_5^1-\lambda T_6^1=\lambda T_4^0+2\left(1-\lambda \right)T_5^0+\lambda T_6^0$

$-0.0835T_4^1+2.167T_5^1-0.0835T_6^1=0$ …… (6)

$-\lambda T_5^1+2\left(1+\lambda \right)T_6^1-\lambda T_7^1=\lambda T_5^0+2\left(1-\lambda \right)T_6^0+\lambda T_7^0$

$-0.0835T_5^1+2.167T_6^1-0.0835T_7^1=0$ …… (7)

$-\lambda T_6^1+2\left(1+\lambda \right)T_7^1-\lambda T_8^1=\lambda T_6^0+2\left(1-\lambda \right)T_7^0+\lambda T_8^0$

$-0.0835T_6^1+2.167T_7^1-0.0835T_8^1=0$ …… (8)

$-\lambda T_7^1+2\left(1+\lambda \right)T_8^1-\lambda T_9^1=\lambda T_7^0+2\left(1-\lambda \right)T_8^0+\lambda T_9^0$

$-0.0835T_7^1+2.167T_8^1-0.0835T_9^1=0$ …… (9)

$\begin{array}{c}-\lambda T_8^1+2\left(1+\lambda \right)T_9^1-\lambda T_{10}^1=\lambda T_8^0+2\left(1-\lambda \right)T_9^0+\lambda T_{10}^0\\ -\lambda T_8^1+2\left(1+\lambda \right)T_9^1-\lambda T_{10}^1=\lambda T_8^0+2\left(1-\lambda \right)T_9^0+\lambda T_{10}^0\\ -0.0835T_8^1+2.167T_9^1-\left(0.0835\times 50\right)=\left(0.0835\times 50\right)\end{array}$

$-0.0835T_8^1+2.167T_9^1=8.35$…… (10)

Further, the system of equations (2), (3),(4), (5), (6), (7), (8), and (9) can be written in matrix form as,

$\left[\begin{array}{ccccccccc}2.167& -0.0835& & & & & & & \\ -0.0835& 2.167& -0.0835& & & & & & \\ & -0.0835& 2.167& -0.0835& & & & & \\ & & -0.0835& 2.167& -0.0835& & & & \\ & & & -0.0835& 2.167& -0.0835& & & \\ & & & & -0.0835& 2.167& -0.0835& & \\ & & & & & -0.0835& 2.167& -0.0835& \\ & & & & & & -0.0835& 2.167& -0.0835\\ & & & & & & & -0.0835& 2.167\end{array}\right]\left[\begin{array}{c}{T}_1^1\\ T_2^1\\ T_3^1\\ T_4^1\\ T_5^1\\ T_6^1\\ T_7^1\\ T_8^1\\ T_9^1\end{array}\right]=\left[\begin{array}{c}16.7\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 8.35\end{array}\right]$

Thus,

$\left[\begin{array}{c}{T}_1^1\\ T_2^1\\ T_3^1\\ T_4^1\\ T_5^1\\ T_6^1\\ T_7^1\\ T_8^1\\ T_9^1\end{array}\right]={\left[\begin{array}{ccccccccc}2.167& -0.0835& & & & & & & \\ -0.0835& 2.167& -0.0835& & & & & & \\ & -0.0835& 2.167& -0.0835& & & & & \\ & & -0.0835& 2.167& -0.0835& & & & \\ & & & -0.0835& 2.167& -0.0835& & & \\ & & & & -0.0835& 2.167& -0.0835& & \\ & & & & & -0.0835& 2.167& -0.0835& \\ & & & & & & -0.0835& 2.167& -0.0835\\ & & & & & & & -0.0835& 2.167\end{array}\right]}^{-1}\times \left[\begin{array}{c}16.7\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 8.35\end{array}\right]$

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,

$\left[\begin{array}{c}{T}_1^1\\ T_2^1\\ T_3^1\\ T_4^1\\ T_5^1\\ T_6^1\\ T_7^1\\ T_8^1\\ T_9^1\end{array}\right]=\left[\begin{array}{c}7.717983\\ 0.2978\\ 0.0114\\ 0.00044\\ 0.00002\\ 0.0002\\ 0.0057\\ 0.1489\\ 3.8589\end{array}\right]$

Furthermore, for the second step to get the values at $t=0.2\text{ s}$, the values of $T_1^1,T_2^1,T_3^1,T_4^1,T_5^1,T_6^1,T_7^1,T_8^1,\text{ and }{T}_9^1$ 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,

$\left[\begin{array}{ccccccccc}2.167& -0.0835& & & & & & & \\ -0.0835& 2.167& -0.0835& & & & & & \\ & -0.0835& 2.167& -0.0835& & & & & \\ & & -0.0835& 2.167& -0.0835& & & & \\ & & & -0.0835& 2.167& -0.0835& & & \\ & & & & -0.0835& 2.167& -0.0835& & \\ & & & & & -0.0835& 2.167& -0.0835& \\ & & & & & & -0.0835& 2.167& -0.0835\\ & & & & & & & -0.0835& 2.167\end{array}\right]\left[\begin{array}{c}{T}_1^2\\ T_2^2\\ T_3^2\\ T_4^2\\ T_5^2\\ T_6^2\\ T_7^2\\ T_8^2\\ T_9^2\end{array}\right]=4\left[\begin{array}{c}7.717983\\ 0.2978\\ 0.0114\\ 0.00044\\ 0.00002\\ 0.0002\\ 0.0057\\ 0.1489\\ 3.8589\end{array}\right]$

Thus, solve as,

$\left[\begin{array}{c}{T}_1^2\\ T_2^2\\ T_3^2\\ T_4^2\\ T_5^2\\ T_6^2\\ T_7^2\\ T_8^2\\ T_9^2\end{array}\right]={\left[\begin{array}{ccccccccc}2.167& -0.0835& & & & & & & \\ -0.0835& 2.167& -0.0835& & & & & & \\ & -0.0835& 2.167& -0.0835& & & & & \\ & & -0.0835& 2.167& -0.0835& & & & \\ & & & -0.0835& 2.167& -0.0835& & & \\ & & & & -0.0835& 2.167& -0.0835& & \\ & & & & & -0.0835& 2.167& -0.0835& \\ & & & & & & -0.0835& 2.167& -0.0835\\ & & & & & & & -0.0835& 2.167\end{array}\right]}^{-1}\times \left[\begin{array}{c}30.8719\\ 1.1913\\ 0.04597\\ 0.001775\\ 0.00010\\ 0.00089\\ 0.0230\\ 0.5956\\ 15.436\end{array}\right]$

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,

$\left[\begin{array}{c}{T}_1^2\\ T_2^2\\ T_3^2\\ T_4^2\\ T_5^2\\ T_6^2\\ T_7^2\\ T_8^2\\ T_9^2\end{array}\right]=\left[\begin{array}{c}14.2889\\ 1.1028\\ 0.0638\\ 0.0033\\ 0.00023\\ 0.00165\\ 0.0319\\ 0.5514\\ 7.1445\end{array}\right]$