The Mathematica code is here.
Consider the heat conduction problem with Neumann (constant flux) at both boundaries of a solid slab. A constant radiant heat flux is imposed on one surface (derivative = -1) and the other surface is thermally insulated (derivative = 0).
If you increase the flux of the heat at one boundary, the temperature inside the system should raise faster over time. However, with the current code, if you increase the flux from 1 to 100, the system behavior doesn't change at all.
In the code, Quiet[NDSolve[...]] is used to suppress the error message: "NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent." The simulation might be incorrect due to this error message.
I want to have the ability to change the system behavior by changing the value of the flux at the boundary.
This is a very important question to me and I couldn't figure it out. Any help will be highly appreciated.
Best Answer
This is really a question about Mathematica, but the short answer is that the code that you're looking at is not correctly handling the boundary and initial conditions.
The PDE in question is the 1-d heat equation on the unit interval, $$\theta_t = \theta_{xx}$$ with Neumann boundary conditions $\theta_x(0,t) = -1$ and $\theta_x(1,t) = 0$ and with initial condition $\theta(x,0) = e^{-1000 x}/1000$.
The problem is that the default method that Mathematica is using to numerically solve this problem does not like the fact that the boundary conditions are inconsistent with the initial condition in the sense that the initial temperature distribution does not satisfy the Neumann condition at $x=1$. The problem is compounded if one changes the boundary condition at $x=0$, as then the initial temperature satisfies neither boundary condition.
This issue is addressed here.