Instability in Method-of-Lines discretization of a 1-D diffusion problem

Dear all,

I'm modeling the emission of chemicals from a building material placed indoor. The model considers the chemical release is driven by the chemical’s internal diffusion within the building material, which is described by the Fick’s 2nd Law. Here’s a diagram and equations of the modeled system:

Deng et al. have published an analytical solution of this model. I will call this “analytical solution” below. https://www.sciencedirect.com/science/article/pii/S1352231003010082

I’m trying to discretize the system in space using method-of-lines to a system of ODEs and then solve it numerically. This method has been described by Yan et al. (Yan 2009 paper) which is called “state-space model”. I will call this “Yan’s method” below. https://www.sciencedirect.com/science/article/pii/S0360132308000711 Here's a diagram of the discretization:

I’m testing many different combinations of chemicals, building materials and indoor configurations using Yan’s method, and compare the results to the analytical solution. For the cases presented in Yan 2009 paper, Yan’s method works well in approximating the analytical solution. However, for certain cases, Yan’s method produces very different results compared to the analytical. Here are two examples:

(1) Formaldehyde in cellulose material with Dm = 4.15e-15 m^2/s (diffusion coefficient), Km = 0.827 (material-air partition coefficient), dm = 0.03 m (material thickness), in a typical North American house. I used two methods to solve the system of ODEs resulting from the method-of-lines discretization: matlab’s ode15s function (Yan’s method-ode15s), and a semi-analytical solution written by my colleague (Yan’s method-lsym). The simulation duration is 100 days. The ode15s chose the time steps automatically, while the lsym used 1 min as time step. The air concentrations and mass fraction emitted from material to air obtained from Yan’s method differ significantly from the analytical solution (see below figure). For the figure, the left two plots are normal scale of y-axis, while the right two plots are log-scale of y-axis. I’ve tried varying the number of layers of discretization from 30 to 200, but it didn’t help.

(2) A hypothetical chemical with Dm = 5.85e-6 m^2/s (diffusion coefficient), Km = 1.86e12 (material-air partition coefficient), dm = 0.01 m (material thickness), in a typical North American house. If I used 199 layers of discretization (N = 200 in attached matlab code line 37), both Yan’s method-ode15s and Yan’s method-lsym are generating mass inside the building material, i.e., the chemical concentrations inside the material become higher than the initial concentration, so the mass fraction emitted becomes negative (see below figure; the blue line in the left bottom plot is also negative but cannot be seen clearly). This is certainly wrong since the building material is the only chemical source in this system. However, it I used only 3 layers of discretization, Yan’s method-lsym generates correct results, although Yan’s method-ode15s still produces negative mass fraction emitted (see below figure).

Thus, it seems that the method-of-lines discretization of my modeled system has some kind of instability related to properties of chemical, building material and indoor configurations: for some cases it works, but for certain cases it doesn’t. However, I cannot figure out what causes this instability. Does anyone have an idea of the source of this instability? Is it due to errors in my matlab codes, or is it real physical instability of the system? I’ve attached my matlab codes here. The two main scripts are “Compare_SS_even_Analyt_NA_house_example_1.m” and “Compare_SS_even_Analyt_NA_house_example_2.m”. Other .m files are functions. Any hints will be appreciated!

[Susi Reddy]

Hi Lei,

I am trying to solve a similar problem, ODE45 works for me but takes a lot of time to solve as the solver chooses small time steps to converge, Did you comeup with any solution to address your issue.

Thanks in advance.