How to implement boundary condition in this case?

Question

Let's start off by NUMERICALLY solving a 1-D steady-state heat transport problem using IMPLICIT FDM.

$DT_{xx}=0; ~T(x=0)=T_{BL}; ~T(x=l)=T_{BR}$

where $D$ is diffusion; $T$ is temperature; subscript $ _{x}$ is the gradient with respective to spatial coordinate $x$.

To implicitly solving the equation, we discretise as follows:

$\begin{align}D\frac{\frac{T_{x+1}-T_{x}}{\Delta x}-\frac{T_{x}-T_{x-1}}{\Delta x}}{\Delta x}=0\\ \frac{D}{\Delta x^{2}}T_{x+1}-\frac{2D}{\Delta x^{2}}T_{x}+\frac{D}{\Delta x^{2}}T_{x-1}=0\end{align}$

At both boundaries, we enforce the boundary condition by allocating a ghost node outside of the boundary. Therefore, on the left

$\begin{align}D\frac{\frac{T_{2}-T_{1}}{\Delta x}-\frac{T_{1}-T_{BL}}{r_{BL}}}{\Delta x}=0\\ -D(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta x r_{BL}})T_{1}+\frac{D}{\Delta x^{2}}T_{2}=\frac{D}{\Delta xr_{BL}}T_{BL}\end{align}$

Where $T_{BL}$ is the temperature applied at the left ghost node (which is a known value), $r_{BL}$ is the VIRTUAL DISTANCE between the first node to the ghost node.

Similarly the right boundary condition can be implemented using $T_{BR}$ (known value);$r_{BR}=r_{BL}$

This equation can be written as $Ax=b$ form solved by tridiagonal solver. Apparently $A$ and $b$ are function of $r_{BL}=r_{BR}$ which one has to fiddle with.

The question come to that : HOW TO IMPLIMENT A PHYSICALLY MEANINGFULL $r_{BL}=r_{BR}$, SO THAT THE SOLUTION WILL NOT BE DISTURBED BY ARBITRARILY CHOOSING VALUES

Actually it is easy to obtain an analytical solution for the problem above:

$T=\frac{T_{BR}-T_{BL}}{l}x+T_{BL}$

and the flux can be calculated by:

$q=-DT_{x}=-D\frac{T_{1}-T_{0}}{l}$

The following figure shows the sensitivity analysis of $r_{BL}=r_{BR}$, and analytical solution are also plotted as comparison (we use $D=0.5,\Delta x=1,l=10,T_{BL}=20,T_{BR}=0$ as example):

one can see that when $r_b$ is between 1e-9 and 1e-2, the result is close enough to the analytical solution. However, if $r_b$ is too large, the result will deviate from right solution ($r_b$=1000). In contrast, if $r_b$ is too small($r_b$=1e-15), the temperature profile may be correct, but the boundary flux will be wrongly calculated.

It is worth to point out that finding suitable $r_b$ in this example problem may not be complicated (still between 1e-9 and 1e-2 is a huge range). However, The real problem comes when time-dependent boundaries into play at transient process. At that stage, a physical based $r_b$ seems to be very important as a fixed $r_b$ may not hold for all time-dependent values.

SOME of my thoughts:

1. I have read some of the ghost boundary problems ghost boundary but havn't found any clear answers on this issue.
2. It seems that this boundary condition may be considered as Robin type($T_{B}+\frac{\partial T_{B}}{\partial x}=constants$). However, it seems there is no way to find out the derivative values (although analytical solution here offers one, but most of the problems may not have analytical solutions).
3. I have heard about this may result from matching significant digits of $T_{b}$ and $T_1$ search 'boundary conductances' in this page. in short, the selected $r_B$ should make $T_{b}$ and $T_1$ partially matched (the first half of the digits match, and the rest half of the digits mismatch). Then the mismatched values are responsible for calculating the outflow. Apparently, if $r_b$ is too big, all the digits are matched, then there is no flow going out. but still such fixed $r_b$ value may be venerable for non-linear boundary conditions.

Any ideas? Thanks in advance!

Answer 1

There is absolutely no reason not to choose $r_{BL}=r_{BR}=\Delta x$. You phrase your question in terms of ghost nodes, but this is a concept one uses for Neumann boundary conditions. It is not necessary for Dirichlet conditions: simply put a node at each $i\Delta x$, $i=0\ldots N$ and in the two equations where the values $T_0, T_N$ appear, you will replace their respective values by the Dirichlet values. That's all that's necessary.

In fact, since your solution is linear, your finite difference scheme must exactly reproduce the exact solution. This is true for a uniform mesh ($r_{BL}=r_{BR}=\Delta x$) as well as any non-uniform choice ($r_{B*}\neq \Delta x$). The fact that you get a solution that is different from what you should points to the fact that there must be a bug in your code.