Software Publications For Students Links About Author
PL
EN
For Students → PMiKNoM Seminars → Class 10
The content of this page is available under the Creative Commons Attribution-ShareAlike License. Valid XHTML 1.0 Transitional Valid CSS! Last modified: 21st of May 2020
PMiKNoM Seminars
Class 1
Class 2
Class 3
Class 4
Class 5
Class 6
Class 7 - Colloquium I
Class 8
Class 9
Class 10
Class 11 and 12
Class 13 - Colloquium II
Class 14 - The retake of the Colloquium
Class 15 - Assessment Colloquium
PMiKNoM Project
Sustainable Development
MSc and BSc theses
Print content

Diffusion-reaction model

PMiKNoM seminars - class 10

Introduction

Model we've discussed on seminar 8 assumed that the transport of chloride ions occurs in the uniform, continuous system. However, in reality the concrete sample is a porous systems, and diffusion of ions occurs primarily in the porous liquid, inside the pores of concrete. In order to include the information about porosity in the model, the left side of Fick's Law must be multiplied by the porosity coefficient `\phi`:

$$ \phi \frac{\partial c}{\partial t} = D \frac{ \partial^2 c}{\partial x^2} $$

Additionally, the adsorption reaction occurs and some of the ions are bound in the solid phase:

$$ \phi \frac{\partial c}{\partial t} = D \frac{ \partial^2 c}{\partial x^2} - r $$
Reaction term `r` represents ions, that are adsorbed in time `\partial t`, and the change in bound chloride concentration `{\partial C_b}/{\partial t}` is described as:
$$ (1 -\phi) \rho_s \frac{\partial C_b}{\partial t} = r $$

In the equation above the density of the solid phase `rho_s` occurs, because `c` and `C_b` are expressed in different units, i.e. [`kg`/`m^3`] and [`kg`/`kg`] respectively.

Adsorption (reaction term) can be described using three isotherms: linear, Langmuir or Freunlich. We will use the latter. Details of the model including equations, numerical solution and programs algorithm can be found below.

Time-dependent diffusion-reaction model - basic equations

Diffusion and binding

After including the information about porisity of concrete and chloride binding in solid phase, the mass ballance equation takes the form:

$$ \left\{ \begin{array}{l} \phi \frac{\partial c}{\partial t} = D \frac{ \partial^2 c}{\partial x^2} - r \\ (1 -\phi) \rho_s \frac{\partial C_b}{\partial t} = r \end{array} \right. $$
where `c` represents the concentration of free chlorides in porous liquid, `C_b` represents the concentration of bound chlorides, `\phi` represents porosity, `D` is the diffusion coeffisient, `x` is space, and `t` represents time.

The reaction term, `r` is calculated using Freundlich isotherm:

$$ r = k [c - (C_b \text{/} K_b)^{1\text{/}\eta}] $$
where `k`, `K_b` and `\eta` are reaction constants.

The density of the solid phase in concrete, `\rho_s` is calculated using the weighted sum:

$$ \rho_s = \frac{\rho_c - \phi \cdot \rho_w}{1 - \phi} $$
where `rho_c` oznaczarepresents total density of concrete, and `rho_w = 1000 {kg}/m^3` is the density of water.

Initial and boundary conditions

We assume the constant concentrations at the boundaries (Dirichlet boundary conditions):

$$ \left\{ \begin{array}{l} c(0,t) = c_L = const \\ c(l,t) = c_R = const \end{array} \right. $$

At the beginning of the process, there is no chlorides in the sample:

$$ \left\{ \begin{array}{l} c(x,0) = 0 \\ C_b(x,0) = 0 \end{array} \right. $$
The total concentration of chlorides

The total concentration of chlorides (rescaled to `kg`/`kg` of concrete) is a sum of free and bound chlorides:

$$ C_{t o t}(x,t) = C_{f r e e}(x,t) + C_{b o u n d}(x,t), $$
where     `C_{f r e e}(x,t) = \frac{\phi \cdot c(x,t)}{\rho_c}`     and     `C_{b o u n d}(x,t) = \frac{(1- \phi) \cdot \rho_s \cdot C_b(x,t)}{\rho_c}`.

Diffusion-reaction model - Numerical solution

Numerical solution of the above equations is based on the Method of Lines (MoL) and on the explicit Euler method.

Method of Lines

We divide the considered segment into `n+1` sub-segments, each having the length `h`, to obtain the discrete points `x_0` to `x_{n+1}`. Let's then define `c_i \equiv c(x_i)` and `C_{b,i} \equiv C_b(x_i)`. The second derivative can be now substituted with the finite difference. That leads to the following set of equations:

$$ \left\{ \begin{array}{lll} \phi \frac{dc_i}{dt} = \frac{D}{h^2} (c_{i+1}-2c_i+c_{i-1}) - k [c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta}] & & \text{ for } i=1..n \\ (1 -\phi) \rho_s \frac{dC_{b,i}}{dt} = k [c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta}] & & \text{ for } i=0..n+1 \end{array} \right. $$
Eulers Method

If `c_{i}` and `C_{b,i}` denote the concentrations at time `t_j`, while `\overline c_{i}` and `\overline C_{b,i}` denote the chloride concentrations at time `t_{j+1}` then:

$$ \left\{ \begin{array}{ll} \overline c_{i} = c_{i} + dt \cdot \frac{dc_{i}}{dt} \\ \overline C_{b,i} = C_{b,i} + dt \cdot \frac{dC_{b,i}}{dt} \end{array} \right. $$

After inserting the time derivative, one obtains:

$$ \left\{ \begin{array}{ll} \overline c_{i} = c_{i} + \alpha \left( c_{i+1}-2c_i+c_{i-1} \right) - \beta \bigl( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} \bigr) \\ \overline C_{b,i} = C_{b,i} + \gamma \bigl( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} \bigr) \end{array} \right. $$
where `\alpha = \frac{dt \cdot D}{\phi \cdot h^2}`, `\beta = \frac{dt \cdot k}{\phi}` and `\gamma = \frac{dt \cdot k}{(1 -\phi) \rho_s}`.

Diffusion-reaction model - programs algorithm

  1. Declare the global constants `c_L`, `c_R`, `l`, `D`, `t_{end}`, `n`, `m`, as well as `row`, `roc`, `fi`, `k`, `Kb` and `eta`.
  2. Begin the main procedure with declarations of variables `h`, `dt`, `ros`, `\alpha`, `\beta`, `\gamma`, `i`, `j`, and arrays `c(0 \text{ to } n+1)`, `c_\text{new}(1 \text{ to } n)` and `C_b(0 \text{ to } n+1)`.
  3. Calculate the step value `h`, the time-step value `dt`, the density `ros` and the variables `\alpha`, `\beta`, `\gamma`.
  4. Write the initial condition for every point: `c(i)=0` and `C_b(i)=0`.
  5. Write the boundary condition: `c(0)=c_L` and `c(n+1)=c_R`.
  6. Execute the loop, in which for every time-step (`j=1` to `m`) three inner loops will be executed:
    • Program calculates the new values of concentration of free chlorides, and assign them to the array `c_\text{new}`, i.e.:
      for `i = 1` to `n`     `c_{\text{new},i} = c_{i} + \alpha ( c_{i+1}-2c_i+c_{i-1} ) - \beta ( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} )`.
    • Program calculates the new values of concentration of bound chlorides, and assign them to the array `C_b`, i.e.:
      for `i = 0` to `n+1`     `C_{b,i} = C_{b,i} + \gamma ( c_i - (C_{b,i} \text{/} K_b)^{1\text{/}\eta} )`.
    • Program switches the arrays, i.e.:
      for `i = 1` to `n`     `c_i=c_{\text{new},i}`.
  7. Display results, in columns containing: `i \cdot h`, `C_{f r e e}(i)`, `C_{b o u n d}(i)`, `C_\{t o t}(i)`. (See equation 2.3.)

Diffusion-Reaction model - the obtained software