An efficient numerical model for liquid water uptake in porous material and its parameter estimation
Ainagul Jumabekova Julien Berger Denys Dutykh Hervé Le Meur Aurélie Foucquier Mickaël Pailha Christophe Ménézo


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

A. Jumabekova]Ainagul Jumabekova
A. Jumabekova: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, UMR 5271 CNRS, LOCIE, 73000 Chambéry, France
Ainagul.Jumabekova@univ-smb.fr
https://www.researchgate.net/profile/Ainagul_Jumabekova/

J. Berger]Julien Berger
J. Berger: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, UMR 5271 CNRS, LOCIE, 73000 Chambéry, France
Berger.Julien@univ-smb.fr
https://www.researchgate.net/profile/Julien_Berger3/

D. Dutykh]Denys Dutykh
D. Dutykh: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAMA, 73000 Chambéry, France and LAMA, UMR 5127 CNRS, Université Savoie Mont Blanc, Campus Scientifique, F-73376 Le Bourget-du-Lac Cedex, France
Denys.Dutykh@univ-smb.fr
http://www.denys-dutykh.com/

H. Le Meur]Hervé Le Meur
H. Le Meur: Laboratoire de Mathématiques d'Orsay, Univ. Paris-Sud, CNRS, Université Paris-Saclay, 91405 Orsay, France
Herve.LeMeur@math.u-psud.fr
https://www.math.u-psud.fr/` lemeur/

A. Foucquier]Aurélie Foucquier
A. Foucquier: Univ. Grenoble Alpes, CEA, LITEN, DTS, INES, F-38000, Grenoble, France
Aurelie.Foucquier@gmail.com
https://www.researchgate.net/profile/Aurelie_Foucquier/

M. Pailha]Mickael Pailha
M. Pailha: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, UMR 5271 CNRS, LOCIE, 73000 Chambéry, France
Mickael.Pailha@univ-smb.fr

Ch. Ménézo]Christophe Ménézo
Ch. Ménézo: Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, UMR 5271 CNRS, LOCIE, 73000 Chambéry, France
Christophe.Menezo@univ-smb.fr
https://www.researchgate.net/profile/Christophe_Menezo/




empty 
Ainagul Jumabekova

LOCIE–CNRS, Université Savoie Mont Blanc, France


Julien Berger

LOCIE–CNRS, Université Savoie Mont Blanc, France


Denys Dutykh

LAMA–CNRS, Université Savoie Mont Blanc, France


Hervé Le Meur

Université Paris-Saclay, Orsay, France


Aurélie Foucquier

CEA–LITEN, DTS, France


Mickael Pailha

LOCIE–CNRS, Université Savoie Mont Blanc, France


Christophe Ménézo

LOCIE–CNRS, Université Savoie Mont Blanc, France




  [t]1.0
    

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

    

  





 





empty 

   

 
Last modified: December 30, 2023





empty





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.



: water uptake process; porous material;  numerical scheme; parameter estimation problem; advection–diffusion equation with gravity flux 


MSC: [2010] 35R30 (primary), 35K05, 80A20, 65M32 (secondary)


PACS: [2010] 44.05.+e (primary), 44.10.+i, 02.60.Cb, 02.70.Bf (secondary)







§ 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 <cit.>. 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 <cit.>.

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 <cit.>. Some examples of such deterioration can be found in <cit.> or <cit.>.

Thus, one understands the importance of studying the physical phenomena occurring in rising damp. As mentioned by Franzoni in <cit.>, 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) <cit.>. 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 <cit.>, 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 <cit.>. However, the main drawback is that these solutions work only for idealized conditions. As an example, in <cit.> 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 <cit.>. 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 <cit.>, 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  numerical scheme. This approach is particularly efficient for so-called advection–diffusion equations as highlighted from a mathematical point of view in <cit.> and illustrated in <cit.> 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 <ref>, the physical model is described. Following Section <ref> presents the  scheme. Its accuracy and advantages are compared to classical methods on numerical case studies. Section <ref> describes the experimental facility used to generate observation of water uptake. Section <ref> 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 <cit.>. 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:

    θt -  ÷ 𝐣 ,

where θ [ 𝗆^ 3/𝗆^ 3 ] is the volume basis of liquid water content in the material and 
 || 𝐣 ||  [ 𝗄𝗀/(𝗆^ 2.𝗌) ] is the liquid flow through the capillaries.

The liquid flow is driven by diffusion, advection due to air pressure difference and gravity forces <cit.>:

    jj_ dj_ aj_ g .


The diffusion liquid transfer flux is expressed as: 

    j_ d D_ θ θ D_ T  T  ,

where D_ θ [ 𝗆^ 2/𝗌 ] is the liquid transport coefficient under a water content gradient and D_ T [ 𝗆^ 2/(𝗌.𝖪) ] is the liquid transport coefficient associated to a temperature gradient. Both depend on water content θ. Since the experiments are conducted under isothermal conditions, the liquid diffusion under temperature gradient is assumed negligible. Thus, we have:

    j_ d D_ θ θ .

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

    j_ g K  i ,

where K  [ 𝗆/𝗌 ] is the hydraulic conductivity depending on θ, and i is taken as pointing upward as illustrated in Figure <ref>. 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:

    j_ aΠ θ  ,

where [ 𝗆/𝗌 ] is the mass average velocity and Π 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: 

    ( 1H/L )  _ 0 ,

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

    H ( t ) ∫_ 0^ L θ ( x,t ) /θ_ sat dx  ,

where  θ_ sat – saturation moisture content.

Let us define coefficient a_ 0 and approximate it using Darcy's law as:

    a_ 0Π _ 0 ≈ Π k_ a/μ Δ P/L ,

where k_ a [ 𝗆^ 2 ] is the intrinsic permeability of the material, and μ [ 𝖯𝖺· 𝗌 ] is the dynamic viscosity of the fluid.

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

    θtx (  D_ θ θx a_ 0  ( 1 H )  θ i K  i )  .


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:

    θ ( x0 ,t )  0  ,    t 0 
    θ_ sat ,    t >  0   , 
            θ ( x L ,t )  0  ,
            θ ( x , t 0)  0  .








 §.§ 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:

    u θ/θ_ sat , 
             x^ ⋆x/L ,
             t^ ⋆t/t_ ref ,

where θ_ sat – saturation moisture content, L – length of the brick, t_ ref –  set to one hour.

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

    d^ ⋆ ( u ) D_ θ( θ )/D_ ref ,
             a_ 0^ ⋆a_ 0/a_ ref ,
                k^ ⋆ ( u ) K( θ )/k_ ref .


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

    Fot_ ref D_ ref/L^ 2 .

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

    Pea_ ref L/D_ ref .

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

    Bok_ ref L/θ_ sat D_ ref .

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

    ut^ ⋆   Fo x^ ⋆(  d^ ⋆ ( u )  ux^ ⋆Pe  a_ 0^ ⋆ ( 1  H^ ⋆ )   u Bo  k^ ⋆ ( u )  )   ,
 
where the dimensionless water front height is given by:

    H^ ⋆∫_ 0^ 1  u ( x^ ⋆ ,t^ ⋆ ) d x^ ⋆ .

Further in article, for the sake of clarity the upper script ^ ⋆ is omitted, and Eq. (<ref>) transforms into:

    ut    · x ( d ( u )  · ux ·   a ( u ) ·  u  · k ( u )  )   ,
 
with the following boundary and initial conditions: 

    u( x0 , t)  0  ,    t 0 
    1 ,    t >  0   , 
             u( x1 , t)  0  ,
             u( x , t 0)  0  .

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  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:

    ut    =  x (  d · ux a · u  k (u)  ) ,         t  >  0  ,   x  ∈ [ 0  ,  1 ] ,

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

    u ( 0,t )  u_ L( t )  ,
              u ( 1,t )  u_ R( t ) .


A brief demonstration of the uniqueness of the solution of Eq. (<ref>) is provided for the linearized equation in Appendix <ref>. It can be remarked that the uniqueness of advection–diffusion equation has been demonstrated in <cit.> 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 <cit.>, the problem considered is linear. In <cit.>, 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 <cit.>.

Let us discretize uniformly the space and time intervals, with the parameters Δ x and Δ t, respectively. The discrete values of function u (x, t) are defined as u_j^n  ≡ u ( x_j, t_n ), where j  ∈ { 1, … ,N } and n ∈ { 1, … ,N_t }.





 §.§ The  numerical scheme


The flux J is denoted by the following expression:

    J  d · ux a · u  k (u)  .

For the cell 𝒞 x  ∈ [ x_ j- ,  x_ j+ ] illustrated in Figure <ref>, the semi–discretization of the equation (<ref>) results as:

    u_ j/ t1/Δ x ( J_ j+1/2^ n  J_ j-1/2^ n ) .

Within the  approach, the assumption is that the flux J is constant on the dual cell 𝒞^ ⋆ [ x_ j , x_ j+1 ]. The following boundary value problem can be written:


    J_ j+1/2^ n    =  d · ux a · u  k^n_j+1/2 , 
            ∀ x  ∈  [ x_ j,  x_ j+1 ] ,
    
      u     =  u^ n_ j ,
             x =  x_ j ,
    
      u     =  u^ n_ j+1 , 
             x =  x_ j+1 .


where k^ n_ j+1/2 is approximated by:

    k^ n_ j+1/2 k  ( 1/2 (  u^ n_ j u^ n_ j+1 ) )  .

It is important to note that problem (<ref>) has two unknowns u ( x , t ) and J_ j+1/2^ n with two constrains at x  x_ j and x  x_ j+1. Therefore, one can obtain the exact solution of Eq. (<ref>) as:

    J_ j+1/2^ n a · (  u^n_j+1 u^n_j  e^a  Δ x/d )/1  e^a  Δ x/d k^n_j+1/2 .

Applying expression (<ref>) to Eq. (<ref>), and using an Euler explicit approach, we obtain the expression to compute u_ j^ n+1:

    u^ n+1_ j u^n_jΔ t/Δ x · [  a · (  u^n_j+1 u^n_j-1  e^a  Δ x/d )/1  e^a  Δ x/d
     
      a · (  1  e^a  Δ x/d )/1  e^a  Δ x/d · u^n_j  k^n_j+1/2 k^n_j-1/2 ]  .

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  function  <cit.>.







 §.§ Properties of  scheme



The main advantage of the  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 <cit.>. The  scheme has first order accuracy over the time and space Ø (Δ x + Δ t) <cit.>. Moreover, the flux is also approximated to the first order Ø ( Δ x ).

For Eq. (<ref>), the Courant–Friedrichs–Lewy (CFL) stability condition for the  numerical scheme combined with an Euler explicit approach, is calculated with the expressions from <cit.>:

    Δ t  · max_  j[ (  a dk/du ( u ( x_ j+ ) ) ) · tanh( (  a dk/du ( u ( x_ j+ ) )  )  · Δ x/2 d)^-1 ]
    ≤ Δ x  .

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

As mentioned before, Eq (<ref>) has two unknowns u ( x , t ) and J_ j+1/2^ n. The solution J_ j+1/2^ n is given by Eq. (<ref>). The exact expression of solution u can be computed for x  ∈ 𝒞^ ⋆: 

    u^ n( x ) u^n_j ^a Δ x/d u^n_j+1/^a  Δ x/d 1u^n_j+1 u^n_j/^a  Δ x/d 1 · ^(a · (x  x_j)d), 
     
      x  ∈ [ x_ j ,  x_ j+1 ]  .






 §.§ Extension for the nonlinear case



In case where the coefficients a, d and k of Eq. (<ref>) all depend on u , we approximate them by using the frozen coefficients assumption on the dual cell [ x_ j ,  x_ j+1 ].

Thus, the solution of the boundary-value problem (<ref>) is written as:

    J_ j+1/2^ n a^ n_ j+1/2 · ( u^ n_ j+1 u^ n_ j ^θ^ n_ j+1/2 )/1 ^ θ^ n_ j+1/2 k^ n_ j+1/2 ,

where 

    θ^ n_ j+1/2  a^ n_ j+1/2 · Δ x/d^ n_ j+1/2 ,

and

    a^ n_ j+1/2 a  ( 1/2(  u^n_j + u^n_j+1 )  )  ,
             d^ n_ j+1/2 d  ( 1/2(  u^n_j + u^n_j+1 )  )  ,
             k^ n_ j+1/2 k  ( 1/2(  u^n_j + u^n_j+1 )  )  .

For the nonlinear case the CFL stability condition for the  numerical scheme is is calculated with the expressions from <cit.>:

    Δ t · max_  j  d_ j max_  j[ p_ j+1/2/d_ j+1/2 tanh(p_ j+1/2 · Δ x/2 d_ j+1/2)^-1 ]  ≤ Δ x   ,

where

    p_ j+1/2 a_ j+1/2 u_ j+1/2 · da/du ( u ( x_ j+ ) ) dk/du ( u ( x_ j+ ) )   .






 §.§ Comparing numerical models


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

    ε_ 2 (x)  ≡ √( 1/N_t∑_j=1^N_t( u^ num( x , t_ j ) u^ ref( x,t_ j ) )^2 ) ,

where N_ t is the number of temporal steps.

The local uniform error is defined as the maximum value of ε_ 2 (x):

    ε_ ∞ ≡ max_x  ∈ [ 0, L_ x ] ε_ 2 (x)   .


The significant correct digits (scd) of the solution is calculated according to <cit.>:

    (u)  ≡  - log_10 *u( x,t_ end )   u^ ref( x,t_ end )/u^ ref( x,t_ end )_ ∞ .


The reference function u^ref ( x, t ) is given by a numerical solution of the differential equation based on the Chebyshev polynomial and adaptive spectral methods and obtained using the function  from the  open source package  <cit.>.





 §.§ 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.





  §.§.§ Linear case


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

    d ( u )  d_ 0 0.05  , 
             a ( u )  a_ 0 0.02  , 
            k ( u )  k_ 1  u  ,

with k_ 1 0.5. It should be noted that the dimensionless numbers are set to unity in Eq. (<ref>). The following boundary and initial conditions are set:

    u ( 0 , t )  0.2 ( 1 cos ( π t )^ 2 )  ,
             u ( 1 , t )  0.3 sin( π  t )^ 2 ,
             u ( x , 0 )  0  .

The time domain is defined as t  ∈ [  0 , 3 ]. According to Equation (<ref>), the  scheme CFL condition is given in this case by: 

    Δ t ≤ tanh( ( a_ 0 k_ 1 ) · Δ x 2 d_ 0 )/a_ 0 k_ 1 · Δ x  .

First, the problem is solved by implementing an adaptive in time approach, using  function  with absolute and relative tolerances set to 10^ - 4. The space discretization parameter is set Δ x  10^ - 2. Given the numerical value of the parameter, the CFL condition is Δ t ≤ 10^ - 3.

Figure <ref> shows the variations of the solution u ( x ,t ) as a function of space and time. It demonstrates a satisfactory agreement between the  numerical solution and the reference one. Figure <ref> highlights the ℒ_ 2 error value ε_ 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  solver.









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

Figure <ref> shows the variation of the error ε_ 2 as a function of Δ x with a fixed time discretization parameter Δ t  10^ -4. The CFL condition Δ t  ⩽ 2.5 · 10^ -4 is respected until  Δ x  5 · 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  scheme. Moreover, the  numerical scheme seems to have a second-order accuracy in space 𝒪( Δ x^ 2 ) which is better than the theoretical results mentioned in Section <ref>. The variations of the computational time with the requested level of accuracy is shown in Figure <ref>. For the accuracy ε_ ∞𝒪 ( 10^-3 ) , the finite–difference is ten-time slower than the  approach combined with the  explicit method. It can be noted that the  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 <ref> displays the variation of the significant correct digits with the computational time. In order to reach scdØ( 1.5 ) the finite–difference scheme requires ten times larger CPU time than the  scheme with  explicit approach. The latter is ten times slower than the one with an adaptive approach.

Figure <ref> presents another error study conducted by varying Δ t while the space discretization is fixed to Δ x  10^ -2. For the given parameters, the CFL condition is Δ t  ⩽ 10^ -3. As noticed in Figure <ref>, a bounded solution can be computed until this condition is respected. It can also be noted that the  scheme with  explicit approach has a greater degree of accuracy. Moreover, the scheme is first–order accurate in time 𝒪 ( t ) as mentioned by the theoretical results.





  §.§.§ Nonlinear case


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

    a ( u )  0.02  0.01 u^ 2 , 
             d ( u )  0.05  0.03 u^ 2 , 
             k ( u )  0.05 u^ 2 .

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

    u ( x , 0 )  0  , 
        u ( 0 ,t )  0.8 sin ( π t / 3 )  0.2 sin ( π t / 5 ) , 
        u ( 1 ,t )  0.5 sin ( π t / 4 )  0.3 sin ( π t / 7 )  .

The solution of this problem is computed for the time horizon t  ∈ [  0  ,  15 ]. The  numerical scheme is used with a space discretization step Δ x  0.01 and an adaptive time step with error tolerances set to 10^ - 4.

As shown in Figure <ref>, the profiles of the solution shows a very satisfactory agreement between the  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  scheme with  explicit approach, (ii)  the central finite–differences approach with  explicit and (iii)  the central finite–differences with adapting time step. As noticed in Table <ref>, the  with the adaptive time step scheme is the faster than  with  explicit approach. The computational time is of the same order between the two adaptive in time approaches ( and central finite–differences). However, the  with adaptive time step approach is ten times more accurate than central finite–difference with adaptive time step.



As mentioned in Section <ref>, the numerical model requires the computation of the coefficients d_ j+, a_ j+ and k_ j+ at the interface of the dual cell [ x_ j ,  x_ j+1 ]. According to Eq. (<ref>), 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+, a_ j+ and k_ j+ are calculated by one the following expressions:


    k_ j+    k  ( u_ j u_ j+1/2 )  ,
     k_ j+   1/2 (   k  (  u_ j ) 
     k  (  u_ j+1 )  )  ,
     k_ j+   2   k  (  u_ j ) ·  k  (  u_ j+1 ) /k  (  u_ j )  k  (  u_ j+1 ) ,
     k_ j+   1/3 (   k  (  u_ j ) 
    √(k  (  u_ j )  ·  k  (  u_ j+1 )) k  (  u_ j+1 )  )  ,
     k_ j+   √(k  (  u_ j )  ·  k  (  u_ j+1 )) , 
     k_ j+   ( 1/2 (   k  (  u_ j )^ 3 k  (  u_ j+1 )^ 3 )  )^1/3 ,
     k_ j+   √(1/2 (   k  (  u_ j )^ 2 k  (  u_ j+1 )^ 2 )  ) .



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 <ref> and Figure <ref>. 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+, a_ j+ and k_ j+ is slightly faster for this particular case.







This section examined and analyzed the  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  scheme with an adaptive time step approach. Therefore, in further sections, the solutions are computed with the  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  ×   22  ×  5   𝖼𝗆 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 <ref>, 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  𝗁. The rising damp is observed on one of the faces 11  ×   22  𝖼𝗆. For this, a camera is placed at 1  𝗆 to take a picture of the brick every 3  𝗆𝗂𝗇. 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  𝖯𝖺 is applied thanks to a chamber with a fan. An illustration of the design is shown in Figure <ref>.  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 K.

The height of the water front is taken in the middle of the brick from the pictures of each experiment every 30  𝗆𝗂𝗇. As shown in Figure <ref>, 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 σ_ h 0.5  𝖼𝗆.

For the case when the brick is exposed to ambient air pressure, the experiment takes t_ max 15 𝗁. For the other experiment, with a pressure difference, it lasts t_ max 7 𝗁. Figure <ref> 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 <ref>, 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 <ref>. 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 <cit.> 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  from the  environment, and polynomials with the highest root mean square error 𝐑-𝐬𝐪𝐮𝐚𝐫𝐞 (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 𝗆, t_ref 1 𝗁, D_ ref  10^ -6 𝗆^ 2/𝗌, k_ ref  10^ -10 𝗆/𝗌 and a_ ref 10^ -9 𝗆/𝗌. The saturated volumetric moisture content θ_ sat, introduced in <cit.>, equals  θ_sat 0.3065  𝗆^ 3/𝗆^ 3. It induces the following dimensionless numbers:

    Fo 0.074   ,
            Pe 2.2 · 10^ -4 ,
            Bo 0.7  .

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

    d ( u )  d_ 4  u^ 4  d_ 1  u  d_ 0 ,

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

    a ( u )  a_ 0 ( 1  H  ) ,

where a_ 0 0.7  , and H ∫_ 0^ 1  u ( x ,t ) dx  .  Last, the dimensionless function k ( u ) of the liquid conductivity part:

    k ( u )  k_ 3  u^ 3 ,

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 𝒪 ( 1 ). For this purpose, the normalized derivative  d_ i · d ( u )d_ i  relative to each coefficient is taken and compared with each other. 

    d_ i · d ( u )d_ i{ d_ 4 u^4 ;   d_ 1 u ;   d_ 0 } .

The results are presented in Figure <ref>. One may conclude that coefficient d_ 4 has more influence on the value of function d ( u ), when liquid uptake occurs u 𝒪( 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 <ref> 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 <ref>, Eq.(<ref>) writes as:

    ut   Fo · x ( d ( u ) · uxPe ·   a ( u ) ·  u Bo · k ( u )  )   ,
 
where

    d ( u ) d_ 4 u^ 4 d_ 1 u  d_ 0 ,       a ( u ) a_ 0 ( 1 ∫_ 0^ 1  u ( x ,t )  dx  ) ,          k ( u ) k_ 3  u^ 3 .

Since literature lacks of consistent data <cit.>, 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 ( a_ 0 ,k_ 3 ,d_ 4 ). 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:

    P ≡ { P_ i }{ a_ 0 ,k_ 3 ,d_ 4 } .

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

A parameter P_ i ∈ P is Structurally Globally Identifiable (SGI), if the following condition is satisfied <cit.>:

    ∀ t  ,      u ( P )  u ( P^ ⋆ )  ⇒ P_ iP^ ⋆_ i .


In the case of our model, if u ( x ,t )  ≡ u^ ⋆ ( x ,t ), then  ut ≡ u^ ⋆t and ux ≡ u^ ⋆x, therefore,

    Fo · x (  d( u ) · uxPe · a( u ) · u 
    Bo · k( u ) ) 
    ≡
    Fo · x ( d^ ⋆( u ) · uxPe · a^ ⋆( u ) · u^ ⋆Bo ·  k^ ⋆( u ) )  .


In addition, straightforward replacements in the Eq. (<ref>) shows us the following:

    utFo · (    ( 4 d_ 4 u^ 3 d_ 1 ) ·ux( d_ 4 u^ 4 d_ 1 u  d_ 0 )
      · ^ 2 u x^ 2
       Pe · a_ 0 · ( 1∫_ 0^ 1  u dx ) · ux 3 Bo · k_ 3 · u^ 2 · ux)  ,

which shows that Equality (<ref>) becomes:

    Fo · (    ( 4 d_ 4 u^ 3 d_ 1 )  · ux(  d_ 4  u^ 4 d_ 1 u  d_ 0 )  · ^ 2 ux^ 2
       Pe · a_ 0 · ( 1 ∫_ 0^ 1  u dx ) · ux 3 Bo ·  k_ 3 · u^ 2 · ux) ≡
       Fo · ( ( 4 d^_,⋆_ 4  u^ 3 d_ 1 )  ·ux( d^ ⋆_ 4  u ^ 4 d_ 1 u  d_ 0 ) · ^ 2 ux^ 2
                             Pe · a^⋆_ 0 ( 1∫_ 0^ 1  u dx ) · ux 3 Bo · k^ ⋆_ 3 ·   u^ 2 · ux) ,

and the final expression is

    ( ( d_ 4 d^_,⋆_ 4 ) · ( 4 u^ 3 u^ 4 ^ 2 ux^ 2 ) Pe · ( a_ 0 a^_,⋆_ 0 ) · ( 1 ∫_ 0^ 1  u dx ) · ux
     3 Bo · (  k_ 3 k^_,⋆_ 3 ) · u^ 2 · ux) ≡ 0


Since that { ( 4 u^ 3 u^ 4 ^ 2 ux^ 2 ) , ( 1 ∫_ 0^ 1  u dx ) · ux ,  u^ 2 · ux}  are independent, then d_ 4, k_ 3 and a_ 0 are SGI <cit.>.





 §.§ 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 <cit.>:

    Y_ P_ i ( t )    σ_ p/σ_ H HP_ iσ_ p/σ_ H ∫_ 0^1  X_ P_ idx  ,    X_ P_ iuP_ i ,

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

To compute the sensitivity coefficients, Equation (<ref>) 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  uk_ 3 obtained by solving the following equation:  

    X_ k_ 3tFo · x ( 
          d ( u ) · X_ k_ 3x( d̃ ( u ) · uxPe · a ( u ) Bo · k̃ ( u ) ) · X_ k_ 3
       Bo u^ 3Pe a_ 0 ∫_ 0^ 1 X_ k_ 3 dx ) ,

where

    d̃ ( u )4 d_4 u^3  d_1 ,      k̃ ( u )3 k_3 u^2 .

For the parameter d_ 4, we define X_ d_ 4  ud_ 4 obtained by solving the following equation:  

    X_ d_ 4tFo · x (     d ( u ) · X_ d_ 4x( d̃ ( u ) · uxPe · a ( u ) Bo · k̃ ( u ) ) · X_ d_ 4
        u^ 4 uxPe a_ 0 ∫_ 0^ 1 X_ d_ 4 dx ) .

And last, for the parameter a_ 0, we define X_ a_ 0  ua_ 0 obtained by solving the following equation:  

    X_ a_ 0tFo · x (     d ( u ) · X_ a_ 0x( d̃ ( u ) · uxPe · a ( u ) Bo · k̃ ( u ) ) · X_ a_ 0
       Pe · ( 1∫_ 0^ 1 u( x ,t ) dx  a_ 0 ∫_ 0^ 1 X_ a_ 0 dx ) · u  ) .

The boundary and initial conditions are the same for each coefficient X_ P_ i:

    X_ P_ i ( x 0 ,t ) 0 ,
             X_ P_ i ( x 1 ,t ) 0 ,
             X_ P_ i ( x ,t 0 ) 0 .

The three differential equations are solved with  numerical scheme. For the investigation of the practical identifiability, the a priori parameters values are used. The time values of Y_ k_ 3  hk_ 3 and Y_ d_ 4  hd_ 4 are shown in Figure <ref>. 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 <ref>. 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 <ref>, 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.





  §.§.§ Methodology to solve the parameter estimation problem


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

    J ( d_ 4 , k_ 3 , a_ 0 ):|| H ( d_ 4 , k_ 3 , a_ 0 )  H_ exp ||_ 2 .

The value of H results from the solution of the direct problem (<ref>) for a given set of parameters ( d_ 4 , k_ 3 , a_ 0 ). The value of H_ exp is given by the measurements from the experimental facility and interpolated on the time grid of the numerical scheme.

The cost function J is minimized using function  from the  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:

    d_ 4 ∈ [ 0, 2 ]  , 
             k_ 3 ∈ [ 0, 2 ]  , 
             a_ 0 ∈ [ 0, 2 ]  .


To estimate the quality of the solution of parameter estimation problem, the normalized Fisher matrix is defined according to <cit.>:

    F   [ F_ ij  ],     ∀ ( i, j )  ∈ { 1,...,N_p } ,   F_ ij 1/σ_ H^ 2 ∫_ 0^ t_max Y_ P_ i  Y_ P_ j dt ,

where Y_ P_ i is the sensitivity coefficient of the solution related to the parameter P_ i, σ_ H the measurement uncertainty and N_ p the number of parameters. The matrix F measures the total sensitivity of the system for the measurements to variations of the entire set of parameters P. Under some assumptions detailed in <cit.>, 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 F  is used to assess the estimation uncertainty by computing an error estimator for the parameter P_ i:

    η_ i√(( F^ -1 )_ i i) .

High values of η_ i indicate a possible high error during the parameter identification process.





  §.§.§ Results and discussion


The estimation process is performed with the following discretization parameters Δ x^ ⋆ 0.05, Δ t^ ⋆ 10^ -2, using the  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. (<ref>) 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. (<ref>) is solved 114 times. Each iteration requires approximately 100  𝗌 to solve the direct problem. It results in a overall performance of 3  𝗁. Considering the results from previous section and particularly Table <ref>, the same process using the standard finite–difference scheme and Euler approach would last for about 7.5  𝗁. For the experimental data with pressure, the calculation of solution u of Eq. (<ref>) is computed almost 270 times. The estimation takes about 7.5  𝗁. This value is comparatively small to the 17.5  𝗁 required with the model based on standard approach. The results highlights that the  numerical scheme, applying to parameter estimation problem, save important computational efforts compared to the standard approach with finite–difference scheme.



Figure <ref> 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 <ref>. 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 <ref>. In addition, as noticed in Table <ref>, there is no significant difference between the a priori and estimated values for parameters. It validated the investigations on practical identifiability carried in Section <ref>, 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 <ref> 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 <ref> 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  𝖯𝖺. The propagation of the water front through the brick seems faster in the case of the pressure difference. As noted in Figure <ref>, 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  <ref>, 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 <ref>.



In Figure <ref>, 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 <ref>. 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:

    ( d ( u ) · ux ·   a_ 0 · ( 1 ∫_ 0^ 1  u ( x ,t ) dx ) · u     · k ( u )  )
       Bi · (  u  u_ ∞ ) ,      x  1  ,

where Bi is a dimensionless surface transfer coefficient and u_ ∞ are the ambient conditions in the facility. It implies to estimate this coefficient Bi which is rather difficult experimentally. It can be estimated by adding the coefficient 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:

    ( d ( u ) · ux ·   a_ 0 · ( 1 ∫_ 0^ 1  u ( x ,t ) dx ) · u  · k ( u )  )     0   , 
     
      x     1  ,

To investigate the importance of modifying the definition of the boundary conditions at x  1, additional computation are performed considering the problem Eq. (<ref>) with the estimated parameters ( d_ 4 , k_ 3 ,a_ 0 ) and the boundary condition Eq. (<ref>). The boundary condition at x  0 and the initial condition are unchanged:

    u( x0 , t)     =   0  ,    t 0 
    1 ,    t >  0   , 
             u( x , t 0)   =  0  .

Figure <ref> displays the difference between the solutions computed with Dirichlet and Neumann boundary conditions. The error ε_ 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 <ref> and <ref>. 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:

    H max { x̂:   u ( x̂ ,t )  ≥ û } ,

where û is a chosen value of the amount of liquid. Figure <ref> compares the numerical solution according to the definition of H with Eq. (<ref>) or (<ref>). 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 (<ref>) 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  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 <ref>. The investigations show that the  numerical model is more accurate and faster than the standard approach based on finite–differences and  explicit scheme. The  numerical model has an explicit formulation avoiding costly sub-iterations at each time step. Since its stability condition is nonlinear and scales with Δ t  ≃ Δ x for large space discretization, the  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 <ref>. 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  𝖯𝖺. Then, the comparison between numerical simulations and the experimental observations of the height of the rising front is carried out in Section <ref>. Since the literature lacks consistent data for the material properties <cit.>, 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 <ref> and <ref>, 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  𝖯𝖺 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  numerical scheme enables to save by 50% the computational cost compared to the standard approach.

Even if this work enhances the efficiency of the  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.





§ 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 (<ref>) can be transformed as

    ut    =  x (  d ux u ( a k( u )/u )  ) ,

then regularized into a general advection–diffusion equation:

    ut    =  x (  d uxã u  ) ,
             x  ∈ Ω[ 0  ,  1  ]  ,   t   >   0  ,

where ã ( u )   (  a k ( u )/u). Further in the Section, we assume the linearized representation ã of ã ( u ) around u_ 0:

    ã a k( u_ 0 )u_ 0 ,
 
where u_ 0 is a given function of ( x ,t ).

The following boundary and initial conditions are considered:

    u ( 0 , t ) ^ L ( t )   , 
             u ( 1 , t ) ^ R ( t )   , 
             u ( x , 0 )  0  .

In order to study the uniqueness of the solution, we use the so-called energy method <cit.>. For that purpose, we assume that two solutions u_ 1 ( x , t ) and u_ 2 ( x , t ) satisfy Eqs. (<ref>) and (<ref>). We define w( x , t ) as:

    w ( x , t )    u_ 1 ( x , t )  u_ 2 ( x , t )  ,

solution of the equation:

    wtx (  d  wxã  w  )  ,

with the boundary and initial conditions:

    w ( 1 , t )  0  , 
             w ( 0 , t )  0  ,
             w ( x , 0 )  0  .


If the energy E ( t )   ∫_Ω  w^ 2 ( x , t )  dx is decreasing and E ( 0 )  0, then the solution of Eqs (<ref>) and (<ref>) is unique :

    Et ⩽  0  .

To study the validity of this condition, both sides of Eq. (<ref>) are multiplied by w and integrated over Ω:

    ∫_Ω w  wt dx
    ∫_Ω  w  x (  d  wxã  w  )  dx  ,

which by differentiating under the integral becomes:

    1/2 Et∫_Ω  w  x (  d  wxã  w  )  dx  ,

By performing an integration by parts, we obtain:

    1/2 Et[  w  (  d  wxã  w  )  ]_ x  0^ x  1∫_Ω (  d  wxã  w  )  wx dx  ,

which becomes using the boundary conditions (<ref>):

    1/2 Et∫_Ω  d  ( wx )^ 2∫_Ω ã  w  wx dx  .

The Cauchy–Schwartz inequality writes:

    | ∫_Ω  f   g  |  ⩽ ( ∫_Ω  f^ 2 dx  )^  ( ∫_Ω  g^ 2 dx  )^  ,

for f and g square–integrable realx–value functions. Thus, the inequality (<ref>) yields to:

    1/2 Et  ⩽ ∫_Ω  d  ( wx )^ 2ã ( ∫_Ω  w^ 2dx  )^  ( ∫_Ω ( wx )^ 2dx  )^  .

For u  ∈  W_ 0^ 1 , p, the Poincaré inequality states:

    ||  u  ||_ L^ p ( Ω ) ⩽  C__ Ω ||   u  ||_ L^ p ( Ω ) ,

where C__ Ω is a constant depending on p and Ω only.  Thus, for a one-dimensional problem and for p  2, the inequality (<ref>) becomes:

    1/2 Et  ⩽ ∫_Ω  d  ( wx )^ 2ã  C__ Ω ∫_Ω ( wx )^ 2dx  ,

which can be re-written as:

    1/2 Et  ⩽ (  d ã  C__ Ω )  ∫_Ω ( wx )^ 2 .

Therefore, we have uniqueness of the solution of Eqs (<ref>) and (<ref>) if k ã  C__ Ω ⩾  0. In other terms, the generalized advection coefficient has to be small compared to the diffusion one to have unique solution.




tocsectionReferences
abbrv



