Debate and discussion of VOCs diffusion coefficient derivation process for the Bodalal model

At present, as more and more people realise the importance of environmental protection, they begin to pay attention on the harm of indoor environmental pollution to the human body [1,2,3,4]. Therefore, many experts and scholars [5,6,7] pay more attention on the emission characteristics of VOCs [8,9,10,11]; among them, the mathematical model established by Bodalal et al. [12] is frequently cited; the diffusion coefficient and distribution coefficient of VOCs in building materials are determined by the double-chamber method. In the model, the diffusion coefficient is used to describe the adsorption rate of VOCs in building materials. The model has usability, but there are some logical errors in the mathematical derivation of the model, which can easily cause trouble to the scholars with weak mathematical foundation. The paper points out the error of model derivation and the cause of the error; the correct derivation processes are given in the literature [13,14,15], and a simple method for calculating diffusion coefficient and distribution coefficient is given in the literature [17,18]. It can help scholars to apply the model in a better way.

Bodalal et al. [12] established a mathematical model to describe VOCs migration in building materials by using the static double-compartment method. Absorption and diffusion of VOCs can be divided into three stages: (i) volatile stage, (ii) diffusion and degradation stage and (iii) stable stage. At the beginning of the experiment, VOCs in chamber 1 were mainly volatile, promoted by a concentration gradient; VOCs diffused through the template to chamber 2, and the concentration of chamber 2 began to rise and reached a certain value that could be ignored. This is the volatile phase. In the middle of the experiment, VOCs in chamber 1 were mainly diffused, as VOCs were absorbed by the sample and self-degraded; the concentration of chamber 1 was constantly decreasing, and VOCs spread through the sample to chamber 2, where the concentration kept rising. This is the stage of diffusion and degradation. After a long period of diffusion, the concentrations of the two chambers reached equilibrium and tended to be equal. This is the stable phase. So, VOCs diffusion process in the experiment can be described by the following partial differential equation: (1) ∂C∂t=D∂2C∂x2 {{\partial C} \over {\partial t}} = D{{{\partial ^2}C} \over {\partial {x^2}}}

Boundary conditions are as follows:

At x = 0, C = k_ec₁

Initial conditions are as follows:

At t = 0, 0 ≤ x ≤ l, C = 0

Based on the experimental analysis, the concentration changes of VOCs in the two laboratories can be described by the following equations: c1=DAV1∫0t(∂C∂x)x=0dtc2=c0−DAV2∫0t(∂C∂x)x=ldt \matrix{ {{c_1} = {{DA} \over {{V_1}}}\int_0^t {{\left( {{{\partial C} \over {\partial x}}} \right)}_{x = 0}}dt} \hfill \cr {{c_2} = {c_0} - {{DA} \over {{V_2}}}\int_0^t {{\left( {{{\partial C} \over {\partial x}}} \right)}_{x = l}}dt} \hfill \cr }

Finally, the diffusion coefficient and the partition coefficient can be calculated by the following equations: (2) characteristic equation: tanq=2αqq2−α2 {\rm{characteristic}}\;{\rm{equation}}:\quad \tan q = {{2\alpha q} \over {{q^2} - {\alpha ^2}}} α=keAlVr=−q2Dl2 \matrix{ {\alpha = {{{k_e}Al} \over V}} \cr {r = {{ - {q^2}D} \over {{l^2}}}} \cr } (3) u=ln(4αq2+α(2+α)) u = ln \left( {{{4\alpha } \over {{q^2} + \alpha (2 + \alpha )}}} \right) r and u are, respectively, the slope and intercept of linear functions as follows: ln(c2−c1)c0=ln(4αq2+α(2+α))−(q2Dl2)t \ln {{({c_2} - {c_1})} \over {{c_0}}} = ln \left( {{{4\alpha } \over {{q^2} + \alpha (2 + \alpha )}}} \right) - \left( {{{{q^2}D} \over {{l^2}}}} \right)t

Among them, D is the diffusion coefficient of VOCs, k_e is the distribution coefficient, C is VOCs concentration on the sample material, c₀ is the initial concentration of VOCs in the high concentration, c_i(i = 1,2) is the VOCs concentration in the low and high concentration compartments, x is the linear distance of VOCs diffusion, t is the diffusion time, l is the thickness of the sample materials, A is the area of the sample and V_i (i = 1,2) is the volume of the chamber.

The model is reasonable, but there are some logical errors in the mathematical derivation of solving the problem. Now the error of model derivation and the cause of the error are explained as follows.

First, the Laplace–Carson transform is applied to the Eq. (1) along with boundary conditions. The result is shown in the following equations: S1(q)=c0qα2αqcosq+α2sinq−q2sinqS2(q)=c0q(αcosq−qsinq)2αqcosq+α2sinq−q2sinq \matrix{ {{S_1}(q) = {{{c_0}q\alpha } \over {2\alpha q\cos q + {\alpha ^2}\sin q - {q^2}\sin q}}} \hfill \cr {{S_2}(q) = {{{c_0}q(\alpha \cos q - q\sin q)} \over {2\alpha q\cos q + {\alpha ^2}\sin q - {q^2}\sin q}}} \hfill \cr }

Among them, q=l(pD)0.5 q = l{\left( {{p \over D}} \right)^{0.5}} , α=keAlV \alpha = {{{k_e}Al} \over V} .

The generalised Heaviside theorem is used to solve the inverse transformation of S2(q)−S1(q)c0 {{{S_2}(q) - {S_1}(q)} \over {{c_0}}} ; so, the relationship between c₁ and c₂ is shown below: (4) c2−c1c0=∑i=1∞4αe−Dqi2l2tqi2+α(α+2) {{{c_2} - {c_1}} \over {{c_0}}} = \sum\limits_{i = 1}^\infty {{4\alpha {e^{ - {{Dq_i^2} \over {{l^2}}}t}}} \over {q_i^2 + \alpha (\alpha + 2)}}

In fact, this is a use of the residue method to solve the inverse transformation of S2(q)−S1(q)c0 {{{S_2}(q) - {S_1}(q)} \over {{c_0}}} . The specific process is that, substituting the q=l(pD)0.5 q = l{\left( {{p \over D}} \right)^{0.5}} , p=q2Dl2 p = {{{q^2}D} \over {{l^2}}} into the S2(q)−S1(q)c0⋅eptp {{{S_2}(q) - {S_1}(q)} \over {{c_0}}} \cdot {{{e^{pt}}} \over p} and solving derivative of its inverse with variable p, so inverse transformation is solved through combining characteristic Eq. (2) and the generalized Heaviside Theorem, but the result is like this (5) c2−c1c0=∑i=1∞4αeDqi2l2tqi2+α(α+2) {{{c_2} - {c_1}} \over {{c_0}}} = \sum\limits_{i = 1}^\infty {{4\alpha {e^{{{Dq_i^2} \over {{l^2}}}t}}} \over {q_i^2 + \alpha (\alpha + 2)}} rather than the results of Eq. (4). Apparently, the infinite series on the right side of Eq. (5) is discrete as time increases, when q is the first positive root, the results cannot be concluded as follows: c2−c1c0=4αe−Dq2l2tq2+α(α+2) {{{c_2} - {c_1}} \over {{c_0}}} = {{4\alpha {e^{ - {{D{q^2}} \over {{l^2}}}t}}} \over {{q^2} + \alpha (\alpha + 2)}} in this way, linear function ln(c2−c1)c0=ln(4αq2+α(2+α))−(q2Dl2)t ln {{({c_2} - {c_1})} \over {{c_0}}} = ln \left( {{{4\alpha } \over {{q^2} + \alpha (2 + \alpha )}}} \right) - \left( {{{{q^2}D} \over {{l^2}}}} \right)t t cannot be given, because of it, the latter result will not occur, there is no way to use this result to calculate the diffusion coefficient. Based on the above analysis, it is necessary to correct the errors.

The right derivation process is given as follows:

The original problem {∂C∂t=D∂2C∂x2x=0, C=kec1x=l, C=kec2t=0, 0≤x≤l, C=0 \left\{ {\matrix{ {{{\partial C} \over {\partial t}} = D{{{\partial ^2}C} \over {\partial {x^2}}}} \hfill \cr {x = 0,\quad C = {k_e}{c_1}} \hfill \cr {x = l,\quad C = {k_e}{c_2}} \hfill \cr {t = 0,\quad 0 \le x \le l,\quad C = 0} \hfill \cr } } \right. is transformed into a boundary value problem of ordinary differential equations with parameter p as follows: {Dd2Sdx2=pSS(0,p)=keS1(p)S(l,p)=keS2(p) \left\{ {\matrix{ {D{{{d^2}S} \over {d{x^2}}} = pS} \hfill \cr {S(0,p) = {k_e}{S_1}(p)} \hfill \cr {S(l,p) = {k_e}{S_2}(p)} \hfill \cr } } \right. by using the Laplace–Carson Transformation, the special solution of the problem can be obtained, (6) S1(q′)=c0q′α2αq′coshq′+α2sinhq′+q′2sinhq′ {S_1}({q^\prime}) = {{{c_0}{q^\prime}\alpha } \over {2\alpha {q^\prime}\cosh {q^\prime} + {\alpha ^2}\sinh {q^\prime} + {q^\prime}^2\sinh {q^\prime}}} (7) S2(q′)=c0q′(αcoshq′−q′sinhq′)2αq′coshq′+α2sinhq′+q′2sinhq′ {S_2}({q^\prime}) = {{{c_0}{q^\prime}(\alpha \cosh {q^\prime} - {q^\prime}\sinh {q^\prime})} \over {2\alpha {q^\prime}\cosh {q^\prime} + {\alpha ^2}\sinh {q^\prime} + {q^\prime}^2\sinh {q^\prime}}}

Among them, q′=lpD {q^\prime} = l\sqrt {{p \over D}} , V₁ = V₂ = V, α1=α2=α=keAlV {\alpha _1} = {\alpha _2} = \alpha = {{{k_e}Al} \over V} .

Using sinh(q′) = −isin(iq′),cosh(q′) = cos(iq′) to organise, the Eqs (6) and (7) take the following form: S1(q)=c0qα2αqcosq+α2sinq−q2sinqS2(q)=c0q(αcosq−qsinq)2αqcosq+α2sinq−q2sinq \matrix{ {{S_1}(q) = {{{c_0}q\alpha } \over {2\alpha q\cos q + {\alpha ^2}\sin q - {q^2}\sin q}}} \hfill \cr {{S_2}(q) = {{{c_0}q(\alpha \cos q - q\sin q)} \over {2\alpha q\cos q + {\alpha ^2}\sin q - {q^2}\sin q}}} \hfill \cr } there, q = iq′, i is imaginary unit.

The p=−q2Dl2 p = - {{{q^2}D} \over {{l^2}}} can be obtained from q=il(pD)0.5 q = il{\left( {{p \over D}} \right)^{0.5}} . Substituting the q=ilpD q = il\sqrt {{p \over D}} into the S2(q)−S1(q)c0⋅eptp=qα(cosq−1)−qsinq2αqcosq+α2sinq−q2sinq⋅eptp {{{S_2}(q) - {S_1}(q)} \over {{c_0}}} \cdot {{{e^{pt}}} \over p} = {{q\alpha (\cos q - 1) - q\sin q} \over {2\alpha q\cos q + {\alpha ^2}\sin q - {q^2}\sin q}} \cdot {{{e^{pt}}} \over p} and solving derivative with variable p, then p=−q2Dl2 p = - {{{q^2}D} \over {{l^2}}} is introduced into the above derivative result. Combining characteristic Eq. (2) and generalised Heaviside theorem, the inverse transformation is solved. The result is written as follows: (8) c2−c1c0=∑i=1∞4αe−Dqi2l2tqi2+α(α+2) {{{c_2} - {c_1}} \over {{c_0}}} = \sum\limits_{i = 1}^\infty {{4\alpha {e^{ - {{Dq_i^2} \over {{l^2}}}t}}} \over {q_i^2 + \alpha (\alpha + 2)}} where q_i is all positive roots of Eq. (2).

The series on the right side of Eq. (8) converges rapidly. So for larger t, an approximate solution of Eq. (8) can be derived by taking the first root of the Eq. (2). That is, c2−c1c0=4αe−Dq2l2tq2+α(α+2) {{{c_2} - {c_1}} \over {{c_0}}} = {{4\alpha {e^{ - {{D{q^2}} \over {{l^2}}}t}}} \over {{q^2} + \alpha (\alpha + 2)}} A linear function can be derived by taking the logarithm of both sides of this equation, as follows: ln(c2−c1)c0=ln(4αq2+α(2+α))−(q2Dl2)t \ln {{({c_2} - {c_1})} \over {{c_0}}} = ln \left( {{{4\alpha } \over {{q^2} + \alpha (2 + \alpha )}}} \right) - \left( {{{{q^2}D} \over {{l^2}}}} \right)t The data of c2−c1c0 {{{c_2} - {c_1}} \over {{c_0}}} and t are determined by experiments; the value of u=ln(4αq2+α(2+α)) u = ln \left( {{{4\alpha } \over {{q^2} + \alpha (2 + \alpha )}}} \right) and r=−q2Dl2 r = {{ - {q^2}D} \over {{l^2}}} are fitted by the least-square regression method. In this way, the diffusion coefficient D=−rl2q2 D = - {{r{l^2}} \over {{q^2}}} and the partition coefficient ke=αVAl {k_e} = {{\alpha V} \over {Al}} can be calculated; among them, the values of r, l, A, V are known, and the values of α and q are calculated by using Eqs (2) and (3). The final result is consistent with the Bodalal model.

The correct derivation is given above. The diffusion and distribution coefficients can be calculated by the following procedure.

First, for the data of c₀, c₁, c₂ from the experiment, the least square method is used to fit the linear relationship between ln lnΔcc0 \ln {{\Delta c} \over {{c_0}}} and t (as shown in Figure 1). Thus, the data of r and u are derived, r is the slope of the linear relationship and u is the intercept of the linear relationship.

Fig. 1

Second, give the scope of α, q according to the reality, find the values of α, q that minimise the function f (α,q) = [α + α cos q − q sin q]² and g(α,q) = {e^u[q² + α(2 + α)] − 4α}² to zero in this range.

In the end, the diffusion coefficient and the distribution coefficient are calculated from the relationships: D=−l2rq2 D = {{ - {l^2}r} \over {{q^2}}} and ke=αVAl {k_e} = {{\alpha V} \over {Al}} .

In the second step, a small programme can be made to solve the values of α, q, for example, the small programme in Figure 2.

Fig. 2

A specific application of the model is given below. First, The VOCs concentration curve with time in white cement putty is obtained through experiments (as illustrated in Figure 3). And then, based on the derivation result of Section 4 and combining the method of Section 5, the value of α, q can be accurately calculated according to the experimental data of u. Finally, in Table 1, computational data of the diffusion coefficient and the distribution coefficient are presented; it is assumed that A = V = l = 1.

Fig. 3

The mathematical model established by Bodalal to describe the physical diffusion process of VOCs has been applied by many scholars in the development of building materials for environmental protection. The improved derivation process in this paper is clearer and can help readers to better understand the model. At the same time, the calculation method of the diffusion coefficient and the distribution coefficient is introduced in this paper, which is simple and easy to be realised by a small programme; it has universality.