2人の大学院生による雑記ブログ

# 7-2　Implicit Method and Finite Volume Method[python]

In Part 7-2, we will discretize using the implicit method in the finite volume method.

As an example, the one-dimensional diffusion equation can be written as:

$$\frac{\partial }{\partial t}(\rho C T)=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)$$

We assume $$T=T(t,x), \lambda=\lambda(x)$$

Discretize in cell j.

We normally discretize the left-hand side time-term with first-order accuracy.

$$\frac{\partial }{\partial t}(\rho C T)=\rho C \frac{\Delta T^n _i}{\Delta t}$$

Attention! $$\Delta T^n = T^{n+1}-T^n$$

Then on the right side, unlike the past, the next time $$n + 1$$ is used instead of the current time $$n$$.

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)= \frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right) ^{n+1}$$

Then we will make a detailed transformation. First, perform Taylor expansion in the time direction as shown below.

$$T^{n+1}=T^n+\Delta t \frac{\partial T}{\partial t} +O(\Delta t^2)$$

$$\simeq T^n+\Delta t \frac{\partial T}{\partial t}$$

$$\simeq T^n +\Delta t \frac{T^{n+1}-T^n}{\Delta t}$$

$$= T^n +\Delta T^n$$

As before，$$\Delta T^n = T^{n+1}-T^n$$

From the above

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right) ^{n+1}= \frac{\partial}{\partial x} \left( \lambda \frac{\partial}{\partial x}( T^n +\Delta T^n)\right)$$

It looks strange but I will tell you how to get it step by step.

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right) =\frac{1}{\Delta x}(\lambda_{i+1/2}\frac{\partial T}{\partial x}_{x=i+1/2}-\lambda_{i-1/2}\frac{\partial T}{\partial x}_{x=i-1/2})$$

$$=\frac{1}{\Delta x} ( \lambda_{i+1/2} \frac{T_{i+1}-T_i}{\Delta x}-\lambda_{i-1/2} \frac{T_{i}-T_{i-1}}{\Delta x} )$$

It means;

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x} \right) ^{n+1}=\frac{1}{\Delta x^2} (\lambda_{i+1/2} (T_{i+1}^n-T_i^n )-\lambda_{i-1/2} (T_{i}^n -T_{i-1}^n))$$ $$+ \frac{1}{\Delta x^2} ( \lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2} (\Delta T_{i}^n -\Delta T_{i-1}^n))$$

Summarize the above results,

$$\frac{\partial }{\partial t}(\rho C T)=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)$$

$$⇔ \rho C \frac{\Delta T^n_i}{\Delta t} = \frac{1}{\Delta x^2}(\lambda_{i+1/2} (T_{i+1}^n-T_i^n )-\lambda_{i-1/2}(T_{i}^n -T_{i-1}^n))$$ $$+ \frac{1}{\Delta x^2}(\lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2}(\Delta T_{i}^n -\Delta T_{i-1}^n))$$

Now we have a discretization for the one-dimensional case.

However, it isn’t finished.

Since $$\Delta T^n= T^{n+1}- T^{n}$$ is a future value, this is an unknown variable and there are three in the above formula.

Bring the unknown variable to the left side then the right side will be known values.

$$⇔ \rho C \frac{\Delta T^n_i}{\Delta t}- \frac{1}{\Delta x^2}(\lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2}(\Delta T_{i}^n -\Delta T_{i-1}^n))$$ $$=\frac{1}{\Delta x^2}(\lambda_{i+1/2} (T_{i+1}^n-T_i^n )-\lambda_{i-1/2}(T_{i}^n -T_{i-1}^n))$$

$$⇔\frac{\lambda_{i-1/2}}{\Delta x^2}\Delta T_{i-1}^n+( \frac{\rho C}{\Delta t}-\frac{\lambda_{i-1/2}}{\Delta x^2}-\frac{\lambda_{i+1/2}}{\Delta x^2})\Delta T_{i}^n+\frac{\lambda_{i+1/2}}{\Delta x^2}\Delta T_{i+1}^n=R_i^n$$

$$\left[ \begin{array}{ccc} \frac{\lambda_{i-1/2}}{\Delta x^2} & ( \frac{\rho C}{\Delta t}-\frac{\lambda_{i-1/2}}{\Delta x^2}-\frac{\lambda_{i+1/2}}{\Delta x^2}) & \frac{\lambda_{i+1/2}}{\Delta x^2} \end{array}\right] \left[ \begin{array}{c} \Delta T_{i-1}^n \\ \Delta T_{i}^n \\ \Delta T_{i+1}^n \end{array}\right] =R_i^n$$

The above formula is about cell i. In other words, it can be created as many as the number of cells.

For example, if there are 10 spatial cells, there are 10 unknown variables.（$$\Delta T ^n_i,i=1～10$$）

And the number of formulas that can be created is 10.

Since there are 10 unknowns and 10 expressions, you can find the answer by combining them.

A matrix is used to solve the simultaneous equations.

The formulas for the entire space can be summarized using a matrix as follows:

$$\left[\begin{array}{ccccc} \ddots& \ddots & \ddots && \\ & \frac{\lambda_{i-1/2}}{\Delta x^2} & ( \frac{\rho C}{\Delta t}-\frac{\lambda_{i-1/2}}{\Delta x^2}-\frac{\lambda_{i+1/2}}{\Delta x^2}) & \frac{\lambda_{i+1/2}}{\Delta x^2} & \\ && \ddots & \ddots & \ddots \\ \end{array}\right] \left[\begin{array}{c} \vdots \\ \Delta T_{i-1}^n \\ \Delta T_{i}^n \\ \Delta T_{i+1}^n \\ \vdots \end{array}\right] =\left[\begin{array}{c} \vdots \\ \vdots \\ R_i^n \\ \vdots \\ \vdots \end{array}\right]$$

$$⇔\bf A \bf \Delta \bf T = \bf R$$

$$⇔\bf \Delta \bf T = \bf A^{-1} \bf R$$

The top and bottom of this matrix $$\bf A$$ will change slightly due to the boundary conditions, but it will look like this.

If you can find the inverse matrix, you will know $$\Delta T^n$$, and the following calculation will give you the final answer.

$$T^{n+1}=T^{n}+\Delta T^n$$

We will show you how to get this inverse matrix in next article.

In the case of two dimensions, the number of coordinates will increase and it will be difficult. Similarly, if the discretization is performed in the cell $$(i,j)$$ by the implicit method, the result will be as follows.

$$\frac{\partial }{\partial t}(\rho C T)=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)+ \frac{\partial}{\partial y} \left( \lambda \frac{\partial T}{\partial y}\right)$$

$$⇔ \rho C \frac{\Delta T^n_i }{\Delta t}=\frac{1}{\Delta x^2}(\lambda_{i+1/2,j} (T_{i+1,j}^n-T_{i,j}^n )-\lambda_{i-1/2,j}(T_{i,j}^n -T_{i-1,j}^n))$$ $$+ \frac{1}{\Delta x^2}(\lambda_{i+1/2,j} (\Delta T_{i+1,j}^n-\Delta T_{i,j}^n )-\lambda_{i-1/2,j}(\Delta T_{i,j}^n -\Delta T_{i-1,j}^n))$$ $$+\frac{1}{\Delta y^2}(\lambda_{i,j+1/2} (T_{i,j+1}^n-T_{i,j}^n )-\lambda_{i,j-1/2}(T_{i,j}^n -T_{i,j-1}^n))$$ $$+ \frac{1}{\Delta y^2}(\lambda_{i,j+1/2} (\Delta T_{i,j+1}^n-\Delta T_{i,j}^n )-\lambda_{i,j-1/2}(\Delta T_{i,j}^n -\Delta T_{i,j-1}^n))$$

This time, it became complicated because the thermal conductivity $$\lambda$$ was not made constant. We want to tell you the way which you can use generally .

Next time, we will calculate the inverse matrix using the iterative method.

*The code can be found in Part 7-5

by hide 