How to implement boundary condition in this case?

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.

• The question, you link to, is on cell centered finite volume methods (F V M), where ghost cells are used for Dirichlet boundaries. For F D M, ghost nodes are used for Neumann conditions. See @WolfgangBangerth 's answer,
– Jan
Oct 22 '13 at 14:37

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.
• if we put node at each $i\Delta x,i=0...N$, the problem domain has $N+1$ nodes. By assigning $T_0$ and $T_N$ as respectively $T_{BL}$ and $T_{BR}$, The unknown node number is only $N-1$ (or equivalently the dimension of x in $Ax=b$ is $N-1$), but if we use Neumann boundary on both sides where the flux between ghost and boundary node are given The unknown node number is $N+1$ (or equivalently the dimension of x in $Ax=b$ is $N+1$). Does it mean under fixed node numbers in domain, one has to change the dimension of $x$ in $Ax=b$ based on the boundaries? which may be a bit annoying at 2d or 3d? Oct 22 '13 at 23:03
• In principle, the answer is yes. In practice, we usually implement in the following way: We just build a linear system with $N+1$ equations, regardless of the type of boundary conditions. Then, if a node lies at a Dirichlet boundary, we modify the corresponding row: we replace it by all zeros except for the diagonal, which we set to one. The right hand side of this equation will then be set equal to $T_{B*}$, guaranteeing that the corresponding $T_0$ or $T_N$ equals this value. (In practice, we often put a different value on the diagonal and also modify the columnbut that is immaterial here.) Oct 23 '13 at 2:33
• thanks, by your way $T_0$ can naturally equals to$T_{B*}$ without changing the dimension of x in linear system. but still there are two questions remain. (1) how to calculate the flux through the Dirichlet boundaries, especially when $(T_0-T_1)/dx\neq (T_1-T_2)/dx$? (2) in what circumstances do we need to change the diagonal and column values at the line of Dirichlet boundaries as you mentioned? I suspect this may be somehow related to the ghost node or equivalently $r_{B*}$ in this question? Oct 23 '13 at 3:46
• About your first question: Yes, you compute the flux as $(T_0-T_1)/dx$. About the second one: the modification of the diagonal element is done for numerical stability reasons -- we typically choose the diagonal element to be of the same order as the other diagonal elements of the matrix, and then scale the right hand side. You can modify the column by doing one Gauss elimination step with this degree of freedom, resulting in a column that is zero with the exception of the diagonal element, if you want to restore symmetry of the matrix. Oct 24 '13 at 11:43