1-D boundary value problem: How implement mixed boundary conditions using a FD method?

I have been given a convection-diffusion ODE modeling the steady state temperature of a pipe (through which flows a fluid) as

$$-\frac{d}{dz}\left(\kappa \frac{dT}{dz} \right)+v\rho C\frac{dT}{dz}=Q(z)$$.

Here $$T$$ is the (steady state) temperature and it depends on the space variable $$z$$. All other parameters are constants and for the function $$Q$$ holds

$$Q(z)= \left\{\begin{matrix} 0,&\quad 0\leq z

As for the boundary conditions (BCs) we have a Dirichlet condition for the leftmost boundary and a Neumann condition for the rightmost boundary,

$$T(0)=T_{0}\quad \text{and}\quad -\kappa \frac{dT}{dz}(L)=\alpha(v)(T(L)-T_{out}),$$

where $$T_{out}$$ is the temperature outside the pipe to the right of $$z=L$$ and $$\alpha$$ is given by

$$\alpha(v) = \sqrt{\frac{v^{2}\rho^{2}C^{2}}{4}+\alpha_{0}^{2}}-\frac{v\rho C}{2}.$$

Now, I would like to discretize this ODE using a second order finite difference method (FDM), so I propose to use the central difference approximation which approximates the first and second derivatives as

\begin{align*} \frac{du}{dz}(z_{i})&=\frac{u(z_{i+1})-u(z_{i-1})}{2h} + O(h^{2})\\ \frac{d^{2}u}{dz^{2}}(z_{i})&=\frac{u(z_{i+1})-2u(z_{i})+u(z_{i-1})}{h^{2}} + O(h^{2}) \end{align*}

Finally, I would like to discretize the Neumann BC, but here things get tricky. If I use a central difference, I will obtain a point outside the domain (ghost point?). Could I use the condition that the temperature outside of the pipe at $$z=L$$ is $$T_{out}$$, or is this false? And also, how do I implement or force the Neumann condition to be fulfilled? Can I simply replace $$dT/dz$$ in the ODE with the BC when approximating the solution at $$z=L$$? This is what I do right now.

Interestingly enough, my solution seems to blow upp when I increase the number of steps, or mesh points, in the discretization of the $$z$$-domain. Look at these results:

The temperature completely sky rockets! Why is this? (If you would like me to post me code, I am more than willing to do so.)

Best regards,

• I think that type right boundary condition is a Robin BC, instead of a Neumann one. How are you solving your system of equations? – nicoguaro Sep 15 '19 at 23:37
• Yes, that is correct, my fault. In Matlab, I allocate a zero matrix for all the equations I obtain from the central difference approximation. Since Q changes over time, so does the rhs of the equations. Essentially, I solve this system Ax=f, and with Matlabs backslash inverse function I get the solution as x=A\f (that is A^{-1}f). And as I stated already when I implement the Robin condition I simply replace the first derivative evaluated at z=L in the ODE, then I combine this with setting the function value of the point outside the domain to T_{out} (which doesnt feel right...) @nicoguaro – SimpleProgrammer Sep 16 '19 at 8:09

Just use the left-sided derivative at the right boundary. Then your Robin boundary condition at the last grid point $$n$$ is expressed as $$\kappa \left( \frac{1}{2} T_{n-2} −2 T_{n-1} + \frac{3}{2} T_n \right) = \alpha h (T_n - T_{out})$$