# Heat diffusion 

- toc:false
- branch: master
- badges: true
- comments: false
- categories: [mathematics, numerical recipes]
- hide: true

-----
Questions:
- How do I use the Forward-Time Centred-Space method (FTCS) to solve the diffusion equation? 

Objectives:
- Apply the Forward-Time Centred-Space method (FTCS) to solve theat heat diffusion equation

-----

### The diffusion equation is an initial value problem

In the previous section we solved a boundary value problem in the form of Laplace's equation. In this section we will look at an initial value problem, which is a more complex type of PDE. An initial value problem is more complex as we are told the starting conditions and then have to predict future behaviour as a function of time.

The example we will use is the one-dimensional diffusion equation:

\begin{equation}
 \frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial x^2}
\end{equation}

In this case we have a variable $\phi(x,t)$ that depends on position $x$ and time $t$ - so can we not solve it in the same way as finding the $\phi(x,y)$ Laplace's equation, which also had two independent variables?

The problem is that we  only have an *initial* condition in the time dimension - we know the value of $\phi(x,t)$ at $t=0$ but we do not typically know the value of $t$ at a later point. In the spatial dimensions we know the boundary conditions at either end of the grid.

Instead, we will use the <mark>Forward-Time Centred-Space method (FTCS)</mark>.

### There are two steps to the Forward-Time Centred-Space method

#### Step one
Use the finite difference method to express the 1D Laplacian as a set of simulatenous equations:

\begin{equation}
\frac{\partial^2\phi}{\partial x^2} = \frac{\phi(x+a,t)+\phi(x-a,t) - 2\phi(x,t)}{a^2}
\end{equation}


where $a$ is the grid spacing.

Substitute this back into the diffusion equation to give a set of simulataneous ODEs:

\begin{equation}
 \frac{d \phi}{d t} = \frac{D}{a^2}(\phi(x+a,t)+\phi(x-a,t)-2\phi(x,t))
\end{equation}

If there are $N$ grid points then Equation 3 corresponds to a set of $N$ equations. It is an ODE as there is derivative with respect to only one variable - time.

#### Step two
We now have a set of simultaneous ODEs for $\phi(x,t)$. So we can use Euler's method to evolve the system forward in time. Euler's method for solving an ODE of the form $\frac{d\phi}{dt} = f(\phi,t)$ has the general form:

\begin{equation}
\phi(t+h) \simeq \phi(t) + hf(\phi, t).
\end{equation}

Applying this to Equation 3 gives:

\begin{equation}
\phi(x,t+h) = \phi(x,t) + h\frac{D}{a^2}(\phi(x+a,t)+\phi(x-a,t)-2\phi(x,t))
\end{equation}




### The FTCS method can be applied using the Python skills we have developed

Consider a 10cm rod of stainless steel initially at a uniform temperature of 20$^\mathrm{o}$ Celsius. The rod is dipped in a hot water bath at 90$^\mathrm{o}$ Celsius at one end, and held in someone's hand at the other. Assume that the hand is at constant body temperature throughout (27$^\mathrm{o}$ Celsius).

The problem can be represented visually as follows:



Our goal is to calculate the temperature profile of the steel as a function of distance $x$ from the cold side to the hot side, and as a function of time. For simplicity let us assume that the rod is perfectly insulated so that heat only moves horizontally; as a result this problem can be modelled as 1-dimensional. Also assume that neighter the hot water bath or the hand change temperature appreciably.

Thermal conduction is described by the diffusion equation (or heat equation in this context)

\begin{equation}
 \frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial x^2},
\end{equation}

where $D$ is the material dependent thermal diffusivity. For steel $D=4.25\times10^{-6}\mathrm{m}^2\mathrm{s}^{-1}$.

First, let's import the libraries we will be using


```python
import numpy as np
import matplotlib.pyplot as plt
```

Now let's fix some of the constants in the problem


```python
L = 0.1 # width of the rod in metres
D = 4.25e-6 # thermal diffusivity of steel
```

And some of the numerical parameters:


```python
N = 100 # number of divisions in the grid
a = L/N # grid spacing
h = 1e-4 # time step
```

Let's specify the boundary conditions and the initial condition


```python
T_left = 27.0  # temperature fixed on left side of rod
T_right = 90.0 # temperature fixed on right side of rod
T_middle = 20.0 # temperature of rod at the beginning
```

And now create an array $T$ to hold the temperature of the rod and an array $T_{\mathrm{new}}$ to calculate the temperature of the rod after evolving through time. Note that there are $N$ grid divisions but $N+1$ grid points (as we are evaluating at the boundary on each edge).


```python
T = np.empty(N+1, float)  
T_new = np.empty(N+1, float)
```

Apply the boundary conditions and initial condition to our rod:


```python
T[0] = T_left
T[-1] = T_right
T[1:-1] = T_middle
T_new[0] = T_left
T_new[-1] = T_right
```


```python
def evolve(T,T_new):
    c = (h*D)/(a*a)
    T_new[1:N] = T[1:N] + c*(T[2:N+1]+T[0:N-1]-2*T[1:N])
    return T, T_new
    
```

Note that in the function `evolve` we implement Equation 5. We make use of Numpy's element-by-element array manipulation to evaluate the $N-1$ equations in a single step.


```python
epsilon = h/100

times=np.array([0.01,0.1,1,10,100])
t_end = times[-1]+epsilon
t=0.0


while t < t_end:
    
    T_new, T = evolve(T, T_new)
    t+=h
    
    for time in times:
        if abs(t-time)<epsilon:
            plt.plot(T,label="{}s".format(round(t,ndigits=3)))
    
plt.legend()
plt.title("Heat diffusion along a 1D rod")
plt.xlabel("Position (cm)")
plt.ylabel("Temperature (degrees centigrade)")
```

---
Keypoints:

- The diffusion equation is an initial value problem
- There are two steps to the Forward-Time Centred-Space (FTCS) method
- The FTCS method can be applied using the Python skills we have developed

---

Do [the quick-test](https://nu-cem.github.io/CompPhys/2021/08/02/Heat-Diffusion-Qs.html).

Back to [Modelling with Partial Differential Equations](https://nu-cem.github.io/CompPhys/2021/08/02/PDEs.html).

---


```python

```

ed599c310f2dcf34ee541887afe383e9e5910376
ipynb
Jupyter Notebook
_notebooks/2021-08-02-Heat-Diffusion.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
_notebooks/2021-08-02-Heat-Diffusion.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
2021-10-06T08:11:43.000Z
2021-11-24T13:46:08.000Z
_notebooks/2021-08-02-Heat-Diffusion.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
2021-11-03T10:11:28.000Z
2021-11-03T10:11:28.000Z
Qwen/Qwen-72B
1. YES
2. YES