MATLAB: Does the PDE toolbox yield negative values for the concentration, when I solve a diffusion equation

Question

I am using the 2015 version of the PDE toolbox for diffusion modeling. I am modeling a layer that is only a few microns thick which I converted to cm on the y-axis, and the x-axis is the width in cm.

The solve parameters I am using are:

Time: linspace(0,72000,72000)

u(t0): 1e-12

Relative Tolerance: 0.01

Absolute Tolerance: 0.001

and the PDE coefficients are:

C=8e-17

A=0

F=0

D=1

I have also made the geometry in the x and y directions equal so that the problem is square but I still see negative values for the concentration (u) data around the side boundaries of the model where I have specified a Neumann=0 (no flux) as well as a few places within the meshing where the diffused substance doesn't reach by the end of the simulation (after 20 hours, 72000 seconds). For diffusion into the geometry I would expect my values to be between 0.1 and 1E-12 with no values being negative. 

What would be the easiest way to ensure all output values were only positive or zero?

 

Answer 1

The diffusion equation is a parabolic equation. The stability and accuracy of the solution is affected by the size of the time steps. The inaccuracies in the temporal integration lead to unphysical results for the concentration. This can be fixed by reducing the time steps to a smaller value.

I have scaled the original equation such that the temporal integration occurs in units of days. The scaling factor has been absorbed in the diffusion constant C (or Deff). Now we can divide the time range into 20,000 steps except that the size of the step is smaller than the original case, where it was 1 unit.

 

The solution for the issue is to change the diffusion constant so that it is integrating in units of days instead of seconds.

This can be accomplished by doing the following steps in the PDE model:

1. Rescale the solver time range from:

linspace(0,72000,72000) to
linspace(0,1,20000)

2. Modify the "C" parameter of the equation to account for the scaling

8E-17/(24*3600)

This resulted in non-negative values for the concentration.

Other reasons that could explain negative values for the concentration include numerical gradients that might arise due to poor mesh quality in certain regions of the domain.