An efficient numerical model for liquid water uptake in porous material and its parameter estimation

TL;DR

This paper introduces an efficient numerical model using the Scharfetter-Gummel scheme to predict liquid water uptake in porous materials, validated against experiments and including parameter estimation for key coefficients.

Contribution

It presents a novel application of the Scharfetter-Gummel scheme for modeling capillary adsorption with improved accuracy and computational efficiency.

Findings

The model accurately predicts liquid uptake in bricks.

The scheme offers relaxed stability conditions and reduced computational cost.

Parameter estimation effectively calibrates moisture diffusivity and hydraulic conductivity.

Abstract

The goal of this study is to propose an efficient numerical model for the predictions of capillary adsorption phenomena in a porous material. The Scharfetter-Gummel numerical scheme is proposed to solve an advection-diffusion equation with gravity flux. Its advantages such as accuracy, relaxed stability condition, and reduced computational cost are discussed along with the study of linear and nonlinear cases. The reliability of the numerical model is evaluated by comparing the numerical predictions with experimental observations of liquid uptake in bricks. A parameter estimation problem is solved to adjust the uncertain coefficients of moisture diffusivity and hydraulic conductivity.

Tables

Table 1. Table 1. Computational time required to compute the solution of the nonlinear case with the different numerical models.

Numerical model	CPU time $\left[𝗌\right]$	CPU time $[\%]$	${𝜺}_{\infty }$
SG with adaptive time step	$105$	$42$	$1.9\cdot {10}^{-4}$
SG with Euler explicit	$249$	$100$	$5.8\cdot {10}^{-5}$
Central FD with adaptive time step	$101$	$40$	$1.1\cdot {10}^{-3}$
Central FD with Euler explicit	$230$	$92$	$2.1\cdot {10}^{-2}$

Table 2. Table 2. Variation of the error of the solution and of the CPU time of the numerical models according to the interpolation formulas Eq. ( 3.7 ) of the coefficients d ​ ( u ) 𝑑 𝑢 d\,(\,u\,)\, , a ​ ( u ) 𝑎 𝑢 a\,(\,u\,) and k ​ ( u ) 𝑘 𝑢 k\,(\,u\,)\, .

Formula	(a)	(b)	(c)	(d)	(f)	(e)	(g)
${𝜺}_{\infty }$	$2\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$	$1.9\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$	$2\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$	$1.9\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$	$2\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$	$1.9\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$	$1.7\cdot {\mathrm{\hspace{0.17em}10}}^{-3}$
CPU time $\left[𝗌\right]$	21.29	25.47	33.77	25.44	34.52	27.48	25.79
Percentage $[\%]$	$61$	$74$	$98$	$74$	$100$	$80$	$75$

Table 3. Table 3. A priori and estimated values of dimensionless coefficients

Parameters	A priori	No pressure		With pressure
	values	Estimated value	Error ${\eta }_i$	Estimated value	Error ${\eta }_i$
$a_{\mathrm{\hspace{0.17em}0}}^{\star }$	0.7	$0.0052$	$\pm \mathrm{\hspace{0.17em}14}$	$0.55$	$\pm \mathrm{\hspace{0.17em}206}$
$k_{\mathrm{\hspace{0.17em}3}}^{\star }$	0.8	$0.8257$	$\pm \mathrm{\hspace{0.17em}0.01}$	$1.3$	$\pm \mathrm{\hspace{0.17em}0.21}$
$d_{\mathrm{\hspace{0.17em}4}}^{\star }$	0.6	$1.0$	$\pm \mathrm{\hspace{0.17em}0.002}$	$1.72$	$\pm \mathrm{\hspace{0.17em}0.02}$
Residual ${\epsilon }_{\mathrm{\hspace{0.17em}2}}$	1	0.05		0.0087

Equations

$$\frac{\partial \theta }{\partial t}=-\nabla \cdot \mathbf{j},$$

$$\mathbf{j}={\mathbf{j}}_{\mathrm{d}}+{\mathbf{j}}_{\mathrm{a}}+{\mathbf{j}}_{\mathrm{g}}.$$

$${\mathbf{j}}_{\mathrm{d}}=D_{\theta }\nabla \theta +D_T\nabla T,$$

$${\mathbf{j}}_{\mathrm{d}}=D_{\theta }\nabla \theta .$$

$${\mathbf{j}}_{\mathrm{g}}=K\mathbf{i},$$

$${\mathbf{j}}_{\mathrm{a}}=\Pi \theta \mathsf{v},$$

$$\displaystyle\boldsymbol{\mathsf{v}}\ =\ \,\biggl{(}\,1\ -\ \frac{H}{L}\,\biggr{)}\,\boldsymbol{\mathsf{v}}_{\,0}\,,$$

$$H(t)={\int }_0^L\frac{\theta (x,t)}{{\theta }_{\mathrm{sat}}}\mathrm{d}x,$$

$$a_0=\Pi {\mathsf{v}}_0\approx \Pi \frac{k_{\mathrm{a}}}{\mu }\frac{\Delta P}{L},$$

$$\displaystyle\frac{\partial\theta}{\partial t}\ =\ \frac{\partial}{\partial x}\,\biggl{(}\,D_{\,\theta}\,\frac{\partial\theta}{\partial x}\ -\ a_{\,0}\,(\,1\ -\ H\,)\,\theta\,\boldsymbol{i}\ -\ K\,\boldsymbol{i}\,\biggr{)}\,.$$

$$\theta (x=0,t)=\{\begin{array}{ll}0,& t=0\\ {\theta }_{\mathrm{sat}},& t>0\end{array},$$

$$u=\frac{\theta }{{\theta }_{\mathrm{sat}}},$$

$$d^{\star }(u)=\text{nicefrac}{D}_{\theta }(\theta )D_{\mathrm{ref}},$$

$$\mathrm{Fo}=\text{nicefrac}{t}_{\mathrm{ref}}{D}_{\mathrm{ref}}{L}^2.$$

$$\mathrm{Pe}=\text{nicefrac}{a}_{\mathrm{ref}}L{D}_{\mathrm{ref}}.$$

$$\mathrm{Bo}=\text{nicefrac}{k}_{\mathrm{ref}}L{\theta }_{\mathrm{sat}}{D}_{\mathrm{ref}}.$$

$$\frac{\partial u}{\partial t^{\star }}$$

$$H^{\star }={\int }_0^{1}u(x^{\star },t^{\star })\mathrm{d}{x}^{\star }.$$

$$\frac{\partial u}{\partial t}$$

$$u(x=0,t)=\{\begin{array}{ll}0,& t=0\\ 1,& t>0\end{array},$$

$$\displaystyle\frac{\partial u}{\partial t}\$$

$$u(0,t)=u_L(t),$$

$$J=d\cdot \frac{\partial u}{\partial x}-a\cdot u-k(u).$$

$$\frac{\text{d}\,u_{\,j}}{\text{d}\,t}\ =\ \nicefrac{{1}}{{\Delta x}}\;\bigg{(}\;J_{\,j+\frac{1}{2}}^{\,n}\ -\ J_{\,j-\frac{1}{2}}^{\,n}\;\bigg{)}\,.$$

$$\displaystyle J_{\,j+\frac{1}{2}}^{\,n}\$$

$$\displaystyle u\$$

$$\displaystyle k^{\,n}_{\,j+\frac{1}{2}}\ =\ k\,\Biggl{(}\,\frac{1}{2}\,\Big{(}\,u^{\,n}_{\,j}\ +\ u^{\,n}_{\,j+1}\,\Big{)}\,\Biggr{)}\,.$$

$$\displaystyle J_{\,j+\frac{1}{2}}^{\,n}\ =\ \ -\ a\,\cdot\;\nicefrac{{\Bigl{(}\,u^{n}_{j+1}\ -\ u^{n}_{j}\,e^{\nicefrac{{a\,\Delta x}}{{d}}}\,\Bigr{)}}}{{1\ -\ e^{\nicefrac{{a\,\Delta x}}{{d}}}}}\ -\ k^{n}_{j+\frac{1}{2}}\,.$$

$$\displaystyle u^{\,n+1}_{\,j}\ =\ u^{n}_{j}\ +\ \frac{\Delta t}{\Delta x}\,\cdot\,\Bigg{[}\;\ -\ a\,\cdot\,\nicefrac{{\Bigl{(}\,u^{n}_{j+1}\ +\ u^{n}_{j-1}\,e^{\nicefrac{{a\,\Delta x}}{{d}}}\,\Bigr{)}}}{{1\ -\ e^{\nicefrac{{a\,\Delta x}}{{d}}}}}\ +\$$

$$\displaystyle a\,\cdot\,\nicefrac{{\Bigl{(}\,1\ +\ e^{\nicefrac{{a\,\Delta x}}{{d}}}\,\Bigr{)}}}{{1\ -\ e^{\nicefrac{{a\,\Delta x}}{{d}}}}}\,\cdot\,u^{n}_{j}\ -\ k^{n}_{j+\frac{1}{2}}\ +\ k^{n}_{j-\frac{1}{2}}\,\Bigg{]}\,.$$

$$\displaystyle\Delta t\,\cdot\,\max_{\,j}\Biggl{[}\,\biggl{(}\,a\ +\ \frac{\mathrm{d}k}{\mathrm{d}u}\,\Bigl{(}\,u\,\bigl{(}\,x_{\,j+{\textstyle{1\over 2}}}\,\bigr{)}\,\Bigr{)}\,\biggr{)}\,\cdot\,\tanh\Biggl{(}\,\biggl{(}\,a\ +\ \frac{\mathrm{d}k}{\mathrm{d}u}\,\Bigl{(}\,u\,\bigl{(}\,x_{\,j+{\textstyle{1\over 2}}}\,\bigr{)}\,\Bigr{)}\,\biggl{)}\,\cdot\,\nicefrac{{\Delta x}}{{2\,d}}\Biggr{)}^{-1}\,\Biggr{]}$$

Full text

Key words and phrases: water uptake process; porous material; Scharfetter–Gummel numerical scheme; parameter estimation problem; advection–diffusion equation with gravity flux

MSC:

PACS:

Last modified:

1. 1 Introduction

2. 2 Physical model

3. 2.1 Liquid transfer

4. 2.2 Dimensionless formulation

5. 3 Numerical method

6. 3.1 The Scharfetter–Gummel numerical scheme

7. 3.2 Properties of Scharfetter–Gummel scheme

8. 3.3 Extension for the nonlinear case

9. 3.4 Comparing numerical models

10. 3.5 Numerical validation

11. 3.5.1 Linear case

12. 3.5.2 Nonlinear case

13. 4 Experimental facility

14. 4.1 Experimental observations

15. 4.2 Defining the material properties

16. 5 Comparison of the numerical predictions with experimental data

17. 5.1 Structural Identifiability

18. 5.2 Practical Identifiability

19. 5.3 Estimation and comparison with the experimental observations

20. 5.3.1 Methodology to solve the parameter estimation problem

21. 5.3.2 Results and discussion

22. 6 Conclusion

23. A Demonstration of uniqueness of the solution

Introduction

Nowadays efficient energy consumption of buildings is one of the most challenging tasks. In case of the historical buildings, this problem has to deal not only with the energy performance but also with the building protection itself [32, 8]. The study of moisture’s amount in the building walls takes the largest part for energy reduction since it can raise the heat loss through the walls. The increase of moisture quantity may occur due to different reasons as wind–driven rain, rising damp or simply vapour diffusion [4].

Among all the potential sources of moisture, rising damp is a major issue since it is an important source of liquid water. It may strongly impact the building structure by modifying the mechanical behaviour of the wall. Moreover, the water from the ground brings dissolved salts. During the evaporation process, these salts precipitate and crystals appear that deteriorate the material [27, 28]. Some examples of such deterioration can be found in [22] or [26].

Thus, one understands the importance of studying the physical phenomena occurring in rising damp. As mentioned by Franzoni in [10], the main challenges related to the rising damp problem is to find accurate, fast and cheap methods to characterize the moisture content during the water uptake of a material or a wall. Several methods to characterize the amount of moisture already exist and can be divided in two main groups: (i) invasive (traditional like weighting–drying, chromatography) and (ii) non-invasive (dielectric, microwave) [17]. But usually, these approaches are very complex procedures requiring extensive technical support. Moreover, the obtained results do not describe an exact situation of the whole building. Thus, as underlined in [8], it is worth investigating other approaches by proposing numerical models to perform simulations and obtain quantitative results.

Several models were proposed in the literature for the simulation of liquid water uptake in a porous material. A first approach is to propose analytical solutions of the problem as for instance in [13, 18]. However, the main drawback is that these solutions work only for idealized conditions. As an example, in [18] a step–function is used for the diffusivity coefficient, neglecting its nonlinear variation with the fields. This assumption is not valid for the practical case study. Other approaches are based on standard numerical methods in [16, 14, 24, 31]. However, these approaches also have an important drawback. When using explicit Euler approach, the standard stability conditions must be respected implying very small values of discretization time grids due to the high nonlinearity of the problem. Thus these numerical models have a high computational cost. These works promote the use of implicit schemes for their unconditional stability properties. However, they require a large number of sub–iterations during the computation to treat nonlinearities, which also lead to an increase in the computational cost of the numerical models. It is of major importance to propose efficient numerical models to represent the physical phenomena of capillary adsorption in a porous material.

When proposing efficient numerical models, one important point is its reliability to represent the physical phenomena. For this, the numerical predictions need to be confronted to experimental observations. In the case of rising damp, there are several material properties involved in the definition of the physical model. To have an accurate numerical prediction, the moisture diffusivity and the hydraulic conductivity must be known precisely. As mentioned in [14], there is a lack of data in the literature for these material properties. Thus, when comparing the numerical predictions to experimental observation, the material properties may be estimated to calibrate the numerical model. This procedure requires solving the parameter estimation problem. The solution of such type of the inverse problem requires the computation of the so-called direct problem several times. Thus, since current numerical models have high computational cost, the estimation of material properties becomes an obstacle.

Therefore, the objective of this article is twofold. First, an innovative numerical model is proposed for capillary adsorption phenomena in a porous material. It is based on the Scharfetter–Gummel numerical scheme. This approach is particularly efficient for so-called advection–diffusion equations as highlighted from a mathematical point of view in [12] and illustrated in [2, 3] for the case of heat and moisture transfer in building porous materials. In our work, the proposed numerical model is compared to the standard methods in the context of capillary adsorption phenomena. The issue is to prove its reliability, accuracy and smaller computational cost. Then, the numerical model is compared to experimental data to evaluate its accuracy to represent the physical phenomena. To calibrate the model, a parameter estimation problem is solved in a good time to determine the material properties of a brick using observed data of water uptake.

The paper is organized in the following way. In Section 2, the physical model is described. Following Section 3 presents the Scharfetter–Gummel scheme. Its accuracy and advantages are compared to classical methods on numerical case studies. Section 4 describes the experimental facility used to generate observation of water uptake. Section 5 aims at comparing the numerical predictions and the experimental observations, as well as estimating the material properties of the brick to calibrate the model.

Physical model

Liquid transfer

Let us present the physical model of liquid water uptake process in a porous material. The physical model is inspired by the one proposed by Philipp and De Vries in 1957 [25]. Since the problem deals with high moisture content, the exchange between vapor and liquid water is supposed negligible. Moreover, it is assumed that no chemical reaction occurs. Thus, the liquid water conservation equation is given by:

[TABLE]

where \theta\ \mathsf{\bigl{[}\,m^{\,3}/m^{\,3}\,\bigr{]}} is the volume basis of liquid water content in the material and

\Bigl{|}\Bigl{|}\,\overrightarrow{\boldsymbol{\mathbf{j}}}\,\Bigr{|}\Bigr{|}\ \mathsf{\Bigl{[}\,kg/(m^{\,2}.s)\,\Bigr{]}} is the liquid flow through the capillaries.

The liquid flow is driven by diffusion, advection due to air pressure difference and gravity forces [15]:

[TABLE]

The diffusion liquid transfer flux is expressed as:

[TABLE]

where D_{\,\theta}\ \mathsf{\bigl{[}\,m^{\,2}/s\,\bigr{]}} is the liquid transport coefficient under a water content gradient and D_{\,T}\ \mathsf{\bigl{[}\,m^{\,2}/(s.K)\,\bigr{]}} is the liquid transport coefficient associated to a temperature gradient. Both depend on water content $\theta\,$. Since the experiments are conducted under isothermal conditions, the liquid diffusion under temperature gradient is assumed negligible. Thus, we have:

[TABLE]

The transfer is also driven by gravity forces described by the following expression:

[TABLE]

where K\ \mathsf{\bigl{[}\,m/s\,\bigr{]}} is the hydraulic conductivity depending on $\theta\,$, and $\boldsymbol{i}$ is taken as pointing upward as illustrated in Figure 1. To take into account the motion due to the gradient pressure, an advection flux is introduced. It represents liquid water transfer induced by the filtration of dry air:

[TABLE]

where \boldsymbol{\mathsf{v}}\ \mathsf{\bigl{[}\,m/s\,\bigr{]}} is the mass average velocity and $\Pi$ the porosity of the material. Since the filtration of dry air only influences the water front at $x\ =\ H\,$, the mass average velocity is expressed as:

[TABLE]

where $H$ is the water height in the brick, calculated as the total moisture content at the current instant $t\,$:

[TABLE]

where $\theta_{\,\mathrm{sat}}$ – saturation moisture content.

Let us define coefficient $a_{\,0}$ and approximate it using Darcy’s law as:

[TABLE]

where k_{\,\mathrm{a}}\,\mathsf{\bigl{[}\,m^{\,2}\,\bigr{]}} is the intrinsic permeability of the material, and \mu\,\mathsf{\bigl{[}\,Pa\cdot s\,\bigr{]}} is the dynamic viscosity of the fluid.

Finally, the $1-$dimensional liquid transfer equation is expressed as:

[TABLE]

In the physical model, Dirichlet boundary conditions are set since the surface convective transfer coefficients are unknown. At the boundary where liquid uptakes, the water content is set to saturation. At the top of the brick, the water content remains equal to the initial condition where no liquid is present in the brick. Thus, the boundary and initial conditions are defined as:

[TABLE]

Dimensionless formulation

For numerical analysis, it is very useful to work with the unitless formulation of the model. The first reason is the computational accuracy. Scaling the problem to the appropriate order, many terms will have the same order, so the effects of numerical errors are minimized when computing the residual. The second reason is that many useful relationships exist between dimensionless numbers, which show how specific things influence the whole system.

The following dimensionless variables for the moisture content, the time and space domains are defined:

[TABLE]

where $\theta_{\,\mathrm{sat}}$ – saturation moisture content, $L$ – length of the brick, $t_{\,\mathrm{ref}}$ – set to one hour.

Then all the material properties are modified considering a reference value:

[TABLE]

The Fourier number measures the relative importance of the heat and mass transfers through the material:

[TABLE]

The Péclet number demonstrates the relative importance of the advection against the diffusion transfer:

[TABLE]

The importance of the gravity forces with the respect to the diffusion is quantified by the Bond number:

[TABLE]

Finally, the following dimensionless equation of liquid transfer is obtained:

[TABLE]

where the dimensionless water front height is given by:

[TABLE]

Further in article, for the sake of clarity the upper script ${}^{\,\star}$ is omitted, and Eq. (2.2) transforms into:

[TABLE]

with the following boundary and initial conditions:

[TABLE]

Next section will introduce how to solve the retrieved equation effectively.

Numerical method

The issue is now to propose an efficient numerical model to compute an accurate solution with a reduced computational cost. For this, the Scharfetter–Gummel numerical scheme for an advection–diffusion equation with gravity flux is introduced. Its mathematical properties are discussed. For the sake of simplicity, the following differential equation is considered:

[TABLE]

where $\,a\,$, $\,d\,$ are constants and coefficient $\,k\,(u)\,$ depends on $u\,$. Dirichlet boundary conditions are taken as

[TABLE]

A brief demonstration of the uniqueness of the solution of Eq. (3.1) is provided for the linearized equation in Appendix A. It can be remarked that the uniqueness of advection–diffusion equation has been demonstrated in [5] considering a nonlinear advective term and a linear diffusive one with Dirichlet–type boundary conditions. Some demonstrations have been made for multi-dimensional systems of advection–diffusion equations also with Dirichlet–type boundary conditions. In [7], the problem considered is linear. In [21], it is quasi–linear since the diffusion coefficients depend on space and time whereas the advective and source terms are nonlinear. Some proofs have been also produced for a system of diffusion equations with diffusion coefficients depending on the fields and Neumann–type boundary conditions in [1].

Let us discretize uniformly the space and time intervals, with the parameters $\Delta x$ and $\Delta t\,$, respectively. The discrete values of function $u\,(x,\,t)$ are defined as $\,u_{j}^{n}\ \equiv\ u\,(\,x_{j},\,t_{n}\,)\,$, where $\,j\,\in\,\{\,1,\,\ldots\,,N\,\}\,$ and $\,n\,\in\,\{\,1,\,\ldots\,,N_{t}\,\}\,$.

The Scharfetter–Gummel numerical scheme

The flux $J$ is denoted by the following expression:

[TABLE]

For the cell \mathcal{C}\ =\ x\ \in\ \Bigl{[}\,x_{\,j-{\textstyle{1\over 2}}}\,,\,x_{\,j+{\textstyle{1\over 2}}}\,\Bigr{]} illustrated in Figure 2, the semi–discretization of the equation (3.1) results as:

[TABLE]

Within the Scharfetter–Gummel approach, the assumption is that the flux $J$ is constant on the dual cell \mathcal{C}^{\,\star}\ =\ \,\Bigl{[}\,x_{\,j}\,,\,x_{\,j+1}\,\Bigr{]}\,. The following boundary value problem can be written:

[TABLE]

where $k^{\,n}_{\,j+\frac{1}{2}}$ is approximated by:

[TABLE]

It is important to note that problem (3.3) has two unknowns $u\,(\,x\,,\,t\,)$ and $J_{\,j+\frac{1}{2}}^{\,n}$ with two constrains at $x\ =\ x_{\,j}\,$ and $x\ =\ x_{\,j+1}\,$. Therefore, one can obtain the exact solution of Eq. (3.3) as:

[TABLE]

Applying expression (3.4) to Eq. (3.2), and using an Euler explicit approach, we obtain the expression to compute $u_{\,j}^{\,n+1}\,$:

[TABLE]

In this work, several approaches are used for the temporal discretization. As presented above Euler explicit approach can be applied. Additionaly, the solution is computed with adaptive time step the Adams–Bashforth–Moulton algorithm, using Matlab™ function ODE113 [29].

Properties of Scharfetter–Gummel scheme

The main advantage of the Scharfetter–Gummel numerical scheme is its explicit formulation of the solution. Thus, the nonlinearities are handled without costly sub–iterations at each time steps as it is the case with implicit approaches. Another significant point, the scheme is well–balanced as well as asymptotically preserving [12]. The Scharfetter–Gummel scheme has first order accuracy over the time and space $\mathcal{O}\,(\Delta x\ +\ \Delta t)$ [11]. Moreover, the flux is also approximated to the first order $\mathcal{O}\,(\,\Delta x\,)\,$.

For Eq. (3.1), the Courant–Friedrichs–Lewy (CFL) stability condition for the Scharfetter–Gummel numerical scheme combined with an Euler explicit approach, is calculated with the expressions from [12]:

[TABLE]

Through this condition, a nonlinear relation stands between $\Delta t$ and $\Delta x\,$. But if large space discretization $\Delta x$ is taken, time and space grid become proportional to each other $\Delta t\ \simeq\ \Delta x\,$. Thus, the stability condition is relaxed compared to the classical central finite–difference scheme combined with Euler explicit approach $\Delta t\ \simeq\ \Delta x^{\,2}\,$.

As mentioned before, Eq (3.3) has two unknowns $u\,(\,x\,,\,t\,)$ and $J_{\,j+\frac{1}{2}}^{\,n}$. The solution $J_{\,j+\frac{1}{2}}^{\,n}$ is given by Eq. (3.3). The exact expression of solution $u$ can be computed for $x\,\in\,\mathcal{C}^{\,\star}\,$:

[TABLE]

Extension for the nonlinear case

In case where the coefficients $a$, $d$ and $k$ of Eq. (3.1) all depend on $u\,,$ we approximate them by using the frozen coefficients assumption on the dual cell \Bigl{[}\,x_{\,j}\,,\,x_{\,j+1}\,\Bigr{]}\,.

Thus, the solution of the boundary-value problem (3.3) is written as:

[TABLE]

where

[TABLE]

and

[TABLE]

For the nonlinear case the CFL stability condition for the Scharfetter–Gummel numerical scheme is is calculated with the expressions from [12]:

[TABLE]

where

[TABLE]

Comparing numerical models

In order to validate the numerical model, the error between the solution $u^{\,\mathrm{num}}\,(x,\,t)$ and a reference one $u^{\mathrm{ref}}\,(x,\,t)$ is evaluated as a function of $x$ according to the formula:

[TABLE]

where $N_{\,t}$ is the number of temporal steps.

The local uniform error is defined as the maximum value of $\varepsilon_{\,2}\,(x)\,$:

[TABLE]

The significant correct digits (scd) of the solution is calculated according to [30]:

[TABLE]

The reference function u^{\mathrm{ref}}\,\big{(}\,x,\,t\,\big{)} is given by a numerical solution of the differential equation based on the Chebyshev polynomial and adaptive spectral methods and obtained using the function pde23t from the Matlab™ open source package Chebfun [6].

Numerical validation

In this Section, the purpose is to perform numerical computations to validate the proposed model. Thus, only dimensionless quantities are considered with no necessary physical meanings.

3.5.1  Linear case

First, we study a case with constant functions $d\,(\,u\,)$, $a\,(\,u\,)$ and a linear function of $k\,(\,u\,)$:

[TABLE]

with $k_{\,1}\ =\ 0.5\,$. It should be noted that the dimensionless numbers are set to unity in Eq. (2.4). The following boundary and initial conditions are set:

[TABLE]

The time domain is defined as t\ \in\ \bigl{[}\,0\,,\,3\,\bigr{]}\,. According to Equation (3.2), the Scharfetter–Gummel scheme CFL condition is given in this case by:

[TABLE]

First, the problem is solved by implementing an adaptive in time approach, using Matlab™ function ODE113 with absolute and relative tolerances set to $10^{\,-\,4}\,$. The space discretization parameter is set $\Delta x\ =\ 10^{\,-\,2}\,$. Given the numerical value of the parameter, the CFL condition is $\Delta t\leqslant 10^{\,-\,3}\,$.

Figure 3 shows the variations of the solution $u\,(\,x\,,t\,)$ as a function of space and time. It demonstrates a satisfactory agreement between the Scharfetter–Gummel numerical solution and the reference one. Figure 4 highlights the $\mathcal{L}_{\,2}$ error value $\varepsilon_{\,2}\,(\,x\,)$ at the order of $10^{\,-\,4}\,$, which validates the numerical model for this case. This result is consistent with the value of the discretization parameter and the tolerances set in the ODE solver.

An analysis of the numerical error is performed by varying one of the discretization parameters $\Delta x$ and $\Delta t$ values, while the other parameter remains constant. Different numerical schemes are compared for the computation of the solution of problem (3.1). The Scharfetter–Gummel scheme is implemented with (i) the adaptive time step approach using ODE113 solver and a tolerance set to $10^{\,-\,4}$ and (ii) the Euler explicit approach. For comparison, a central finite–difference scheme using the Euler explicit method is also used.

Figure 5(a) shows the variation of the error $\varepsilon_{\,2}$ as a function of $\Delta x$ with a fixed time discretization parameter $\Delta t\ =\ 10^{\,-4}$. The CFL condition $\Delta t\ \leqslant\ 2.5\cdot 10^{\,-4}$ is respected until $\Delta x\ =\ 5\cdot 10^{\,-3}\,$. Beyond this limit, the numerical model cannot provide a bounded solution. These results also confirm a higher accuracy of the numerical model using the Scharfetter–Gummel scheme. Moreover, the Scharfetter–Gummel numerical scheme seems to have a second-order accuracy in space $\mathcal{O}(\,\Delta x^{\,2}\,)\,$ which is better than the theoretical results mentioned in Section 3.2. The variations of the computational time with the requested level of accuracy is shown in Figure 5(b). For the accuracy $\varepsilon_{\,\infty}\ =\ \mathcal{O}\,(\,10^{-3}\,)\,,$ the finite–difference is ten-time slower than the Scharfetter–Gummel approach combined with the Euler explicit method. It can be noted that the Scharfetter–Gummel scheme with the adaptive time step is the fastest one thanks to the nonlinearity of the stability condition and adaptive time step. The computational time can be reduced by a hundred times. Figure 5(c) displays the variation of the significant correct digits with the computational time. In order to reach $\mathrm{scd}\ =\ \mathcal{O}(\,1.5\,)\,$ the finite–difference scheme requires ten times larger CPU time than the Scharfetter–Gummel scheme with Euler explicit approach. The latter is ten times slower than the one with an adaptive approach.

Figure 6 presents another error study conducted by varying $\Delta t$ while the space discretization is fixed to $\Delta x\ =\ 10^{\,-2}\,$. For the given parameters, the CFL condition is $\Delta t\ \leqslant\ 10^{\,-3}\,$. As noticed in Figure 6, a bounded solution can be computed until this condition is respected. It can also be noted that the Scharfetter–Gummel scheme with Euler explicit approach has a greater degree of accuracy. Moreover, the scheme is first–order accurate in time $\mathcal{O}\,(\,t\,)\,$ as mentioned by the theoretical results.

3.5.2  Nonlinear case

The previous validation is performed for a linear case in order to verify the theoretical results and advantages of the Scharfetter–Gummel numerical scheme. Since, the parameter estimation problem to be solved in Section 5 considers nonlinear coefficients, this Section aims at validating the numerical model for such cases. We consider the problem (2.4) with the following nonlinear coefficients:

[TABLE]

The dimensionless numbers are set to unity. The initial and boundary conditions are:

[TABLE]

The solution of this problem is computed for the time horizon t\,\in\,\bigl{[}\,0\,,\,15\,\bigr{]}\,. The Scharfetter–Gummel numerical scheme is used with a space discretization step $\,\Delta\,x\ =\ 0.01\,$ and an adaptive time step with error tolerances set to $10^{\,-\,4}\,$.

As shown in Figure 7, the profiles of the solution shows a very satisfactory agreement between the Scharfetter–Gummel solution and the reference one. It validates the proposed numerical model. The computational time of the numerical model is compared to three other ones: (i) the Scharfetter–Gummel scheme with Euler explicit approach, (ii) the central finite–differences approach with Euler explicit and (iii) the central finite–differences with adapting time step. As noticed in Table 1, the Scharfetter–Gummel with the adaptive time step scheme is the faster than Scharfetter–Gummel with Euler explicit approach. The computational time is of the same order between the two adaptive in time approaches (Scharfetter–Gummel and central finite–differences). However, the Scharfetter–Gummel with adaptive time step approach is ten times more accurate than central finite–difference with adaptive time step.

As mentioned in Section 3.3, the numerical model requires the computation of the coefficients $\,d_{\,j+{\textstyle{1\over 2}}}\,$, $\,a_{\,j+{\textstyle{1\over 2}}}\,$ and $\,k_{\,j+{\textstyle{1\over 2}}}\,$ at the interface of the dual cell \,\bigl{[}\,x_{\,j}\,,\,x_{\,j+1}\,\bigr{]}\,. According to Eq. (3.6), these coefficients are interpolated using the mean values of $\,u_{\,j}\,$ and $\,u_{\,j+1}\,$. The accuracy of this interpolation is tested by using different expressions in the computation of the solution $\,u\,$. The coefficients $\,d_{\,j+{\textstyle{1\over 2}}}\,$, $\,a_{\,j+{\textstyle{1\over 2}}}\,$ and $\,k_{\,j+{\textstyle{1\over 2}}}\,$ are calculated by one the following expressions:

[TABLE]

For each interpolation expression, the CPU time to compute the solution and the error with the reference solution are analyzed. The results are synthesized in Table 2 and Figure 8. All interpolation expressions have a similar order of accuracy. Nevertheless, it highlights that the approach using the arithmetic mean between $u_{\,j}$ and $u_{\,j+1}\,$ to compute the coefficients $\,d_{\,j+{\textstyle{1\over 2}}}\,$, $\,a_{\,j+{\textstyle{1\over 2}}}\,$ and $\,k_{\,j+{\textstyle{1\over 2}}}\,$ is slightly faster for this particular case.

This section examined and analyzed the Scharfetter–Gummel scheme applying to the diffusion–advection equation with linear and nonlinear coefficients. These results demonstrate a significant reduction of the computational cost without losing the accuracy while using the Scharfetter–Gummel scheme with an adaptive time step approach. Therefore, in further sections, the solutions are computed with the Scharfetter–Gummel combined with an adaptive time step approach. Since the numerical model is validated, the next section introduces the experimental facility employed to obtain the observations used in the parameter estimation problem.

Experimental facility

This section presents the facility to produce the experimental observations that will be compared to the numerical predictions of the model to compare its reliability to represent the physical phenomena.

Experimental observations

To perform the liquid water uptake tests, a refractory brick measuring $11\,\times\,22\,\times\,5\ \mathsf{cm}$ is used. The brick is initially sanded allowing a better visualization of the rise of the water. The vertical faces are then taped to be airtight and watertight. As illustrated in Figure 9, a tray is filled with distilled water and maintained at a constant level thanks to a bottle returned with a pierced cap. The brick is placed vertically on a support at the surface of the water. The experiment begins at the moment when the underside of the brick is in contact with the water and lasts at least $7\ \mathsf{h}\,$. The rising damp is observed on one of the faces $11\,\times\,22\ \mathsf{cm}$. For this, a camera is placed at $1\ \mathsf{m}$ to take a picture of the brick every $3\ \mathsf{min}\,$. The experiment occurs in a box opaque to the light to control the luminosity as well as to minimize the convective and radiative exchanges. Some LED lamps provide a constant lighting in the box.

At the top of the brick, two different boundary conditions are imposed. For the first set of experiments, the top of the brick is in contact with the open air at atmospheric pressure. For the second set, a relative pressure of $-\,50\ \mathsf{Pa}$ is applied thanks to a chamber with a fan. An illustration of the design is shown in Figure 10(a). With this device, the air velocity in the brick is assumed as constant in time. The experiment is taken at the constant room temperature $\;T\ =\ 293\,\mathrm{K}$.

The height of the water front is taken in the middle of the brick from the pictures of each experiment every $30\ \mathsf{min}\,$. As shown in Figure 10(b), graph paper is attached to the brick in order to measure the height front. With this setting, the uncertainty measurement of the height scales with $\sigma_{\,h}\ =\ 0.5\ \mathsf{cm}\,$.

For the case when the brick is exposed to ambient air pressure, the experiment takes $\;t_{\,\mathrm{max}}\ =\ 15\;\mathsf{h}\,$. For the other experiment, with a pressure difference, it lasts $\;t_{\,\mathrm{max}}\ =\ 7\;\mathsf{h}\,$. Figure 11 shows the obtained experimental data for both cases.

To verify that the liquid front is homogeneous inside the brick and not occurring only at the interfaces where the pictures are taken, a complementary experiment is performed. The brick is preliminary cut in the middle and then fulfilled with silicone to avoid the liquid uptake in this area. The liquid uptake experiment is performed. At the end of the experiments, the cut enables to rapidly break the brick without perturbing the liquid front. As illustrated in Figure 12, it is confirmed that the liquid front is homogeneous inside the brick. The white silicone that fulfills the cut in the middle of the brick can be noticed.

Defining the material properties

The physical model representing the liquid water uptake test is described in Section 2. Now, it is important to define the material properties involved in the definition of the diffusion, gravity and advection fluxes for the brick used in the experiments. These properties are uncertain [14] and will be determined in the next Section using the experimental data obtained with the facility.

Data from the literature are used to define the dependency of these properties on the moisture content. The data is projected on the polynomial functions, using the cftool from the Matlab™ environment, and polynomials with the highest root mean square error $\mathbf{R-square}$ (the square of the correlation between the response values and the predicted response values) are selected. These measurements provide a priori values of the coefficients of the polynomials, which act as initial guess when solving the parameter estimation problem. The polynomial functions have been chosen for the material properties to fit the data from the literature. Another set of functions can be applied, and the same steps of parameter estimation are utilized. It should be noticed that data fitting is performed only for the liquid water content.

The reference parameters used for computations are $L\ =\ 0.22\;\mathsf{m}\,$, $t_{\mathrm{ref}}\ =\ 1\;\mathsf{h}\,$, $D_{\,\mathrm{ref}}\ =\ 10^{\,-6}\ \mathsf{m^{\,2}/s}\,$, $k_{\,\mathrm{ref}}\ =\ 10^{\,-10}\ \mathsf{m/s}$ and $a_{\,\mathrm{ref}}\ =\ 10^{\,-9}\ \mathsf{m/s}\,$. The saturated volumetric moisture content $\,\theta_{\,\mathrm{sat}}\,$, introduced in [20], equals $\theta_{\mathrm{sat}}\ =\ 0.3065\ \mathsf{m^{\,3}/m^{\,3}}$. It induces the following dimensionless numbers:

[TABLE]

With given reference values, dimensionless representation $d\,(\,u\,)\,$ of the liquid transport coefficient is:

[TABLE]

where $d_{\,4}\ =\ 0.6\,$, $d_{\,1}\ =\ -0.04$ and $d_{\,0}\ =\ 0.0067\,.$ The dimensionless function $a\,(\,u\,)$ of the advection part:

[TABLE]

where $a_{\,0}\ =\ 0.7\,,$ and $H\ =\ \int_{\,0}^{\,1}\,u\,(\,x\,,t\,)\,\mathrm{d}x\,.$ Last, the dimensionless function $k\,(\,u\,)$ of the liquid conductivity part:

[TABLE]

where $k_{\,3}\ =\ 0.8\,.$

Let us discuss which coefficient of $\,d_{\,0}\,$, $d_{\,1}\,$ or $d_{\,4}\,$ contributes the most in the context of liquid uptake or specifically $u\ =\ \mathcal{O}\,(\,1\,)\,$. For this purpose, the normalized derivative $\displaystyle{d_{\,i}\,\cdot\,\frac{\partial d\,(\,u\,)}{\partial d_{\,i}}}$ relative to each coefficient is taken and compared with each other.

[TABLE]

The results are presented in Figure 13. One may conclude that coefficient $d_{\,4}$ has more influence on the value of function $d\,(\,u\,)$, when liquid uptake occurs $u\ =\ \mathcal{O}(\,1\,)\,$. Therefore, in this article only coefficient $d_{\,4}$ is estimated.

Comparison of the numerical predictions with experimental data

The efficiency of the proposed numerical model is demonstrated in Section 3.5 by comparison to reference solution. Now, the purpose is to assess the reliability of the numerical model to predict the physical phenomena. For this, the experimental observations described in previous section are used for comparison with the numerical predictions. With the given parametrization in Section 4.2, Eq.(2.4) writes as:

[TABLE]

where

[TABLE]

Since literature lacks of consistent data [14], the material properties, represented by functions $k\,(\,u\,)\,$, $d\,(\,u\,)$ and $a\,(\,u\,)\,$, are uncertain. Thus, parameter estimation problem are solved to determine the uncertain parameters and calibrate the numerical model with the experimental data. Within the procedure of parameter estimation, first we demonstrate the structural and practical identifiability of the three unknown parameters \bigl{(}\,a_{\,0}\,,k_{\,3}\,,d_{\,4}\,\bigr{)}\,. Then, the parameter estimation problem is solved and the reliability of the numerical model is discussed.

Structural Identifiability

The aim of this section is to demonstrate the formal identifiability of the unknown parameters. The set of unknown parameters is defined as:

[TABLE]

We assume here, $u\,(\,x\,,t\,)$ is the only observable field.

A parameter $\mathrm{P}_{\,i}\,\in\,\mathrm{P}$ is Structurally Globally Identifiable (SGI), if the following condition is satisfied [34]:

[TABLE]

In the case of our model, if $u\,(\,x\,,t\,)\ \equiv\ u^{\,\star}\,(\,x\,,t\,)$, then $\displaystyle\frac{\partial u}{\partial t}\ \equiv\ \frac{\partial u^{\,\star}}{\partial t}$ and $\displaystyle\frac{\partial u}{\partial x}\ \equiv\ \frac{\partial u^{\,\star}}{\partial x}$, therefore,

[TABLE]

In addition, straightforward replacements in the Eq. (5.1) shows us the following:

[TABLE]

which shows that Equality (5.2) becomes:

[TABLE]

and the final expression is

[TABLE]

Since that \displaystyle{\Bigg{\{}\,\Biggl{(}\,4\,u^{\,3}\ +\ u^{\,4}\,\frac{\partial^{\,2}u}{\partial x^{\,2}}\,\Biggr{)}\,,\Biggl{(}\,1\ -\ \int_{\,0}^{\,1}\,u\,\mathrm{d}x\,\Biggr{)}\,\cdot\,\frac{\partial u}{\partial x}\,,\,u^{\,2}\,\cdot\,\frac{\partial u}{\partial x}\Bigg{\}}} are independent, then $\,d_{\,4}\,$, $k_{\,3}$ and $a_{\,0}\,$ are SGI [23].

Practical Identifiability

As demonstrated in previous Section, the three unknown parameters $\,d_{\,4}\,$, $k_{\,3}$ and $a_{\,0}\,$ are globally identifiable. This result remains theoretical. Before solving the parameter estimation problem, it is of major importance to study the sensitivity coefficients of the parameters to demonstrate the practical identifiability. If the sensitivity coefficients are either small or correlated, the estimation problem is difficult and very sensitive to measurement errors. The sensitivity coefficient is defined as the first derivative of the (numerical) observations with respect to an unknown parameter [9, 35]:

[TABLE]

where $\sigma_{\,H}$ and $\sigma_{\,p}$ are scaling factors. In this case, $\sigma_{\,H}$ corresponds to the uncertainty on the experimental observations. The quantity $\sigma_{\,p}$ is set to unity since all unknown parameter have the same uncertainty on the a priori values.

To compute the sensitivity coefficients, Equation (5.1) is differentiated with respect to each unknown parameter. Three differential equations are obtained enabling to compute the sensitivity coefficients. For the parameter $k_{\,3}\,$, we define $X_{\,k_{\,3}}\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\displaystyle\frac{\partial u}{\partial k_{\,3}}$ obtained by solving the following equation:

[TABLE]

where

[TABLE]

For the parameter $d_{\,4}\,$, we define $X_{\,d_{\,4}}\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\displaystyle\frac{\partial u}{\partial d_{\,4}}$ obtained by solving the following equation:

[TABLE]

And last, for the parameter $a_{\,0}\,$, we define $X_{\,a_{\,0}}\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\displaystyle\frac{\partial u}{\partial a_{\,0}}$ obtained by solving the following equation:

[TABLE]

The boundary and initial conditions are the same for each coefficient $\,X_{\,P_{\,i}}\,$:

[TABLE]

The three differential equations are solved with Scharfetter–Gummel numerical scheme. For the investigation of the practical identifiability, the a priori parameters values are used. The time values of $Y_{\,k_{\,3}}\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\displaystyle\frac{\partial h}{\partial k_{\,3}}$ and $Y_{\,d_{\,4}}\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\displaystyle\frac{\partial h}{\partial d_{\,4}}\,$ are shown in Figure 14(a). One may conclude that the estimation of the parameters $k_{\,3}$ and $d_{\,4}\,$ will likely have a good result since their sensitivity coefficients have large magnitudes and they are not correlated with each other. On the other hand, the time variation of the third coefficient $Y_{\,a_{\,0}}$ is given in Figure 14(b). It is not correlated with other coefficients, demonstrating that it is possible to estimate it. However, its magnitudes are very small showing that its estimation with satisfactory accuracy is a difficult task. One can argue that these observations are due to the a priori values for the coefficient $a_{\,0}.$ However, as noticed in Figure 14(b), even multiplying by 15 or 100 times the a priori value of $a_{\,0}\,$, the magnitude of the sensitivity coefficient almost doesn’t change.

Estimation and comparison with the experimental observations

Since the formal and practical identifiability of the three unknown parameters $d_{\,4}\,$, $a_{\,0}$ and $k_{\,3}$ have been demonstrated, the aim is now to estimate them using the experimental observations obtained with the facility. First, the methodology to solve the inverse problem is briefly described. Then, the results are presented and the reliability of the numerical model predictions are discussed.

5.3.1  Methodology to solve the parameter estimation problem

The parameter estimation problem is solved by minimizing the following cost function by the least squares method:

[TABLE]

The value of $H$ results from the solution of the direct problem (5.1) for a given set of parameters \bigl{(}\,d_{\,4}\,,\,k_{\,3}\,,\,a_{\,0}\,\bigr{)}. The value of $H_{\,\mathrm{exp}}$ is given by the measurements from the experimental facility and interpolated on the time grid of the numerical scheme.

The cost function $\mathrm{J}$ is minimized using function fmincon from the Matlab™ environment. This method uses the interior-point algorithm with constraints on the unknown parameters. Here, the upper and lower boundary constraints for the parameters are defined based on the preliminary calculations:

[TABLE]

To estimate the quality of the solution of parameter estimation problem, the normalized Fisher matrix is defined according to [19, 33]:

[TABLE]

where $Y_{\,P_{\,i}}$ is the sensitivity coefficient of the solution related to the parameter $\mathrm{P}_{\,i}$, $\sigma_{\,H}$ the measurement uncertainty and $N_{\,p}$ the number of parameters. The matrix $\mathrm{F}$ measures the total sensitivity of the system for the measurements to variations of the entire set of parameters $\mathrm{P}$. Under some assumptions detailed in [35], the inverse of the Fisher matrix is the matrix of variance of the parameters considered as random variables of the given observable fields. In other words, it summarizes the quality of the information obtained in the parameter identification process. Thus, the inverse matrix of $\mathrm{F}$ is used to assess the estimation uncertainty by computing an error estimator for the parameter $P_{\,i}\,$:

[TABLE]

High values of $\eta_{\,i}$ indicate a possible high error during the parameter identification process.

5.3.2  Results and discussion

The estimation process is performed with the following discretization parameters $\Delta\,x^{\,\star}\ =\ 0.05\,$, $\Delta\,t^{\,\star}\ =\ 10^{\,-2}\,$, using the Scharfetter–Gummel numerical scheme and the Adams–Bashforth–Moulton time adaptive algorithm. The solution of the optimization problem requires the computation of the solution $u$ of Eq. (2.4) for each set of the parameters $(\,d_{\,4},k_{\,3},a_{\,0}\,)$. During the optimization process for the case without pressure, the direct problem Eq. (2.4) is solved $114$ times. Each iteration requires approximately $100\ \mathsf{s}$ to solve the direct problem. It results in a overall performance of $3\ \mathsf{h}\,$. Considering the results from previous section and particularly Table 1, the same process using the standard finite–difference scheme and Euler approach would last for about $7.5\ \mathsf{h}\,$. For the experimental data with pressure, the calculation of solution $u$ of Eq. (2.4) is computed almost $270$ times. The estimation takes about $7.5\ \mathsf{h}\,$. This value is comparatively small to the $17.5\ \mathsf{h}$ required with the model based on standard approach. The results highlights that the Scharfetter–Gummel numerical scheme, applying to parameter estimation problem, save important computational efforts compared to the standard approach with finite–difference scheme.

Figure 15 shows the time variation of the water height for both experiments with pressure and without pressure. It compares the solution computed using the estimated values of the parameters and the experimental data. The results demonstrate a good agreement with the experimental observations, while the solution with the a priori values cannot approximate the experimental data well. One may say that the calibrated numerical model has a satisfactory reliability to represent the physical phenomena of liquid water uptake. Similar conclusions can be noted for the experiments with pressure.

As a result of the optimization process, the estimated values, the residuals of the cost function and the error estimators are reported in the Table 3. The error estimator is small for parameters $k_{\,3}\,$ and $d_{\,4}\,$, proving that the accuracy of estimation is decent. However, the high values of the error estimator of parameter $a_{\,0}\,$ indicate an unsatisfactory estimation with a very high uncertainty. These results are consistent with ones obtained when analyzing the sensitivity coefficient with a priori values in Section 5.2. In addition, as noticed in Table 3, there is no significant difference between the a priori and estimated values for parameters. It validated the investigations on practical identifiability carried in Section 5.2, since it is eligible in the local area of a priori values.

In order to demonstrate the convergence of the optimization algorithm, the following procedure is implemented. The parameter $a_{\,0}\,$ cannot be estimated accurately so its value is fixed for this analysis. Then several values of initial guesses for parameters $d_{\,4}\,$ and $k_{\,3}\,$ are taken. The estimation procedure is performed for each case. Figure 16 shows the path from each starting point until the estimated parameters values. As one may conclude the final points are in a similar range for all the tested initial guesses. It validates the convergence of the optimization method. It should be noted that the figure does not show the intermediary computations of the algorithm procedure.

Figure 17 displays the volume of liquid water relative to space and time when the brick is exposed to the open air or to a pressure difference of $-\,50\ \mathsf{Pa}\,$. The propagation of the water front through the brick seems faster in the case of the pressure difference. As noted in Figure 15, the experimental data with pressure show that the velocity of water uptake is faster since the same height of the water is achieved twice faster than with the experiment without pressure. One may conclude that the pressure influences the velocity of water uptake. However, as reported in Table  3, results from the optimization problem gives different values for the set of parameters $(\,d_{\,4}\,,k_{\,3}\,,a_{\,0}\,)$ with and without pressure. Particularly, values for parameters $\,d_{\,4}\,$ and $\,k_{\,3}\,$ vary from one case to another. It may indicate that material properties may be different between the two bricks used during the experiments. Moreover, the error estimators of the advection coefficients are really poor indicating that for the given value of pressure difference, it has a minor influence on the water uptake. It should also be noted that the sensitivity coefficients of parameter $a_{\,0}$ has really low magnitudes of variations as shown in Figure 14(b).

In Figure 15, it can be noticed that the full length of the brick is not reached, validating the chosen physical model with Dirichlet boundary conditions presented in Section 2. One may argue that other types of boundary conditions could have been chosen at the top of the brick $x\ =\ 1\,$. The Robin-type boundary condition is avoided. Indeed the boundary flow is proportional to the ambient and surface conditions:

[TABLE]

where $\mathrm{Bi}$ is a dimensionless surface transfer coefficient and $u_{\,\infty}$ are the ambient conditions in the facility. It implies to estimate this coefficient $\mathrm{Bi}$ which is rather difficult experimentally. It can be estimated by adding the coefficient $\mathrm{Bi}$ in the inverse problem, at the price of additional computational costs. In addition, the expected values of the surface coefficient are slow and the boundary conditions can be modified into a homogeneous Neumann one:

[TABLE]

To investigate the importance of modifying the definition of the boundary conditions at $x\ =\ 1\,$, additional computation are performed considering the problem Eq. (5.1) with the estimated parameters $(\,d_{\,4}\,,k_{\,3}\,,a_{\,0}\,)$ and the boundary condition Eq. (5.3.2). The boundary condition at $x\ =\ 0$ and the initial condition are unchanged:

[TABLE]

Figure 18 displays the difference between the solutions computed with Dirichlet and Neumann boundary conditions. The error $\varepsilon_{\,2}$ is lower than $10^{\,-4}$ and one may conclude that the boundary conditions do not change significantly the solution of the numerical model. These results are consistent when analyzing Figures 17(c) and 17(a). It can be observed that the flow at $x\ =\ 1$ remains null.

Another open question is the definition of the liquid height. Previously, the height is defined by the integral of $u$ over the space domain. An alternative definition is to the height of the water front in the brick using a threshold:

[TABLE]

where $\hat{u}$ is a chosen value of the amount of liquid. Figure 19 compares the numerical solution according to the definition of $H$ with Eq. (2.3) or (5.4). The numerical predictions are varying with the definition. On one hand, the integral definition is a more general approach but it assumes that the water vapor is negligible in the material. On the other hand, the definition (5.4) requires to determine the value of the threshold. For these reasons, the definition of the liquid height in the brick according remains an open question.

Conclusion

The moisture in the walls greatly affects the overall building performance. Such a process as rising damp has a strong impact on the energy consumption of the building. In this work, the physical model represents the water uptake process in a single brick, based on diffusion, advection and gravity fluxes. To investigate these physical phenomena it is important to have efficient numerical models in terms of computational cost and accuracy. These models require also to have a good reliability to represent the physical phenomena. To answer this issue, an innovative approach is proposed based on the Scharfetter–Gummel numerical scheme. The first part of the paper studies the properties of the numerical model such as accuracy, stability conditions and computational time with a reference solution for both linear and nonlinear cases in Section 3.5. The investigations show that the Scharfetter–Gummel numerical model is more accurate and faster than the standard approach based on finite–differences and Euler explicit scheme. The Scharfetter–Gummel numerical model has an explicit formulation avoiding costly sub-iterations at each time step. Since its stability condition is nonlinear and scales with $\Delta t\ \simeq\ \Delta x$ for large space discretization, the Scharfetter–Gummel numerical scheme combined with an adaptive time step approach is particularly efficient to save computational efforts.

Along with numerical methods, the reliability of the proposed numerical model is analyzed by comparing the numerical predictions to experimental observations. The experimental facility is presented in Section 4. Empirical data are gathered through a set of two experiments. First, the height of the rising front is obtained for a brick under normal conditions. Then, the height front is measured in another brick submitted to a pressure of $-\,50\ \mathsf{Pa}\,$. Then, the comparison between numerical simulations and the experimental observations of the height of the rising front is carried out in Section 5. Since the literature lacks consistent data for the material properties [14], the uncertain parameters are determined by solving a parameter estimation problems. Before solving the parameter estimation problem, the structural and practical identifiability of the three unknown parameters are demonstrated in Sections 5.1 and 5.2, respectively. The results of the parameter estimation problem show that the diffusivity and liquid conductivity are estimated with satisfactory accuracy. The advection coefficient cannot be estimated with accuracy even with a pressure difference of $-\,50\ \mathsf{Pa}$ at the top of the brick. The sensitivity of the numerical solution with respect to this parameter is much smaller compared to the diffusivity and liquid conductivity. With the estimated parameters, there is a satisfactory agreement between the numerical predictions and the experimental observations, highlighting a good reliability of the model. Important computational efforts are saved thanks to the efficiency of the numerical model. The Scharfetter–Gummel numerical scheme enables to save by $50\%$ the computational cost compared to the standard approach.

Even if this work enhances the efficiency of the Scharfetter–Gummel numerical scheme for the solution of advection–diffusion equation with gravity flux, the rising damp problem require simulation of the physical phenomena in $2-$dimensions. Thus, further research should be conducted to extend the numerical model in this way.

Acknowledgements

This work was partly funded by the “Conseil Savoie Mont Blanc” (CSMB), the French Atomic and Alternative Energy Center (CEA) and the French Environment and Energy Management Agency (ADEME) (through the research program CAPVENT). The authors also acknowledge the Junior Chair Research program “Building performance assessment, evaluation and enhancement” from the University of Savoie Mont Blanc in collaboration with the French Atomic and Alternative Energy Center (CEA) and Scientific and Technical Center for Buildings (CSTB). The authors also would like to acknowledge Dr. B. Rysbaiuly for his valuable discussions.

Appendix A Demonstration of uniqueness of the solution

This section proposes to demonstrate the uniqueness of the solution of liquid uptake in porous media. Since $\;u\,(\,x\,,t\,)\,>\,0\;$ at $\;t\,>\,0\,,$ equation (3.1) can be transformed as

[TABLE]

then regularized into a general advection–diffusion equation:

[TABLE]

where \tilde{a}\,(\,u\,)\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\biggl{(}\,a\ +\ \frac{k\,(\,u\,)}{u}\biggr{)}\,. Further in the Section, we assume the linearized representation $\tilde{a}$ of $\tilde{a}\,(\,u\,)\,$ around $u_{\,0}\,$:

[TABLE]

where $u_{\,0}$ is a given function of $\,(\,x\,,t\,)\,$.

The following boundary and initial conditions are considered:

[TABLE]

In order to study the uniqueness of the solution, we use the so-called energy method [7]. For that purpose, we assume that two solutions $u_{\,1}\,(\,x\,,\,t\,)$ and $u_{\,2}\,(\,x\,,\,t\,)$ satisfy Eqs. (A.1) and (A.2). We define $w(\,x\,,\,t\,)$ as:

[TABLE]

solution of the equation:

[TABLE]

with the boundary and initial conditions:

[TABLE]

If the energy $E\,(\,t\,)\,\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\,\displaystyle\int_{\Omega}\,w^{\,2}\,(\,x\,,\,t\,)\,\mathrm{d}x$ is decreasing and $E\,(\,0\,)\ =\ 0\,$, then the solution of Eqs (A.1) and (A.2) is unique :

[TABLE]

To study the validity of this condition, both sides of Eq. (A.3) are multiplied by $w$ and integrated over $\Omega\,$:

[TABLE]

which by differentiating under the integral becomes:

[TABLE]

By performing an integration by parts, we obtain:

[TABLE]

which becomes using the boundary conditions (A.4):

[TABLE]

The Cauchy–Schwartz inequality writes:

[TABLE]

for $f$ and $g$ square–integrable realx–value functions. Thus, the inequality (A.5) yields to:

[TABLE]

For $u\,\in\,W_{\,0}^{\,1\,,\,p}\,$, the Poincaré inequality states:

[TABLE]

where $C_{{}_{\,\Omega}}$ is a constant depending on $p$ and $\Omega\,$ only. Thus, for a one-dimensional problem and for $p\ =\ 2\,$, the inequality (A.6) becomes:

[TABLE]

which can be re-written as:

[TABLE]

Therefore, we have uniqueness of the solution of Eqs (A.1) and (A.2) if $k\ -\ \tilde{a}\,C_{{}_{\,\Omega}}\,\geqslant\,0\,$. In other terms, the generalized advection coefficient has to be small compared to the diffusion one to have unique solution.