Application `feelpp_hm_heat_moisture` devlopment

As is it coded yet, the application feelpp_hm_heat_moisture can only be used for one of the previous benchmarks at a time (drying-layer, and moisture-uptake) : the values of the parameters are distinguished with macros, and we have to build again if we want to change, which is not a good strategy !

The idea is now to adapt the code to make it work for all types of problems we want to model, using a JSON file.

1. Complete variationnal problem

We considere the problem of heat-moisture transport (cf here) :

\[\begin{align} \left(\rho C_{p}\right)_{\mathrm{eff}} \frac{\partial T}{\partial t}+\nabla \cdot\left(-k_{\mathrm{eff}} \nabla T-L_{\mathrm{v}} \delta_{\mathrm{p}} \nabla\left(\phi p_{\mathrm{sat}}\right)\right)=Q \tag{$E_t$}\label{eq:T} \\ \quad \xi \frac{\partial \phi}{\partial t}+\nabla \cdot\left(-\xi D_{\mathrm{w}} \nabla \phi-\delta_{\mathrm{p}} \nabla\left(\phi p_{\mathrm{sat}}\right)\right)=G \tag{$E_\phi$}\label{eq:phi} \end{align}\]

This system of equations is satisfied on the domain $\Omega$. We add those boundary counditions :

\[\begin{align} T &= d_T & \text{on}\,\Gamma_{T,D}\\ \phi &= d_\phi &\text{on}\,\Gamma_{\phi,D}\\ \nabla T\cdot n &= N_T &\text{on}\,\Gamma_{T,N}\\ \nabla\phi\cdot n &= N_\phi &\text{on}\,\Gamma_{\phi,N}\\ \end{align}\]

We have $\partial\Omega=\Gamma_{T,D}\cup\Gamma_{T,N}$, and $\partial\Omega=\Gamma_{\phi,D}\cup\Gamma_{\phi,N}$.

We make the usual time derivate approximation : $\dfrac{\partial T}{\partial t}\approx\dfrac{T^{n+1}-T^n}{\Delta t}$

Let $v$ be a test function associated to $T$. The first equation $\ref{eq:T}$ becomes, after a partial integration :

\[\begin{align} \int_\Omega (\rho C_p)_\mathrm{eff}^n\dfrac{T^{n+1}}{\Delta t}v &+ \int_\Omega(k_\mathrm{eff}^n\nabla T^{n+1})\nabla v + \int_{\Gamma_{T,N}}\left(-k_\mathrm{eff}^n\nabla T^{n+1}\cdot\mathbf{n}\right)v + \int_{\Gamma_{T,D}}\left(-k_\mathrm{eff}^n\nabla T^{n+1}\cdot\mathbf{n}\right)v\\ &+ \int_\Omega L_v^n\delta_\mathrm{p}^np_\mathrm{sat}^n\nabla\phi^{n+1}\nabla v + \int_{\Gamma_{\phi,N}}\left(-L_v^n\delta_\mathrm{p}^np_\mathrm{sat}^n\nabla\phi^{n+1}\cdot\mathbf{n}\right)v + \int_{\Gamma_{\phi,D}}\left(-L_v^n\delta_\mathrm{p}^np_\mathrm{sat}^n\nabla\phi^{n+1}\cdot\mathbf{n}\right)v \\ &+ \int_\Omega L_v^n\delta_\mathrm{p}^n\nabla p_\mathrm{sat}^n \phi^{n+1}\nabla v + \int_{\Gamma_{\phi,N}}\left(-L_v^n\delta_\mathrm{p}^n\nabla p_\mathrm{sat}^n\phi^{n+1}\cdot\mathbf{n}\right)v + \int_{\Gamma_{\phi,D}}\left(-L_v^n\delta_\mathrm{p}^n\nabla p_\mathrm{sat}^n\phi^{n+1}\cdot\mathbf{n}\right)v\\ &=\int_\Omega \left((\rho C_p)_\mathrm{eff}^n\dfrac{T^n}{\Delta t} + Q^{n+1}\right)v \end{align}\]

As $v=0$ on $\Gamma_{T,D}$ and on $\Gamma_{\phi,D}$. There are also Neumann conditions, this equation become :

\[a(T^{n+1},v)+b(\phi^{n+1},v)=f(v)\]

with :

$a(T^{n+1},v) = \displaystyle\int_\Omega (\rho C_p)_\mathrm{eff}^n\dfrac{T^{n+1}}{\Delta t}v + \int_\Omega(k_\mathrm{eff}^n\nabla T^{n+1})\nabla v$
$b(\phi^{n+1},v) = \displaystyle\int_\Omega L_v^n\delta_\mathrm{p}^np_\mathrm{sat}^n\nabla\phi^{n+1}\nabla v + \int_\Omega L_v^n\delta_\mathrm{p}^n\nabla p_\mathrm{sat}^n \phi^{n+1}\nabla v - \int_{\Gamma_{\phi,N}}\left(L_v^n\delta_\mathrm{p}^n\nabla p_\mathrm{sat}^n\phi^{n+1}\cdot\mathbf{n}\right)v$
$f(v) = \displaystyle\int_\Omega \left((\rho C_p)_\mathrm{eff}^n\dfrac{T^n}{\Delta t} + Q^{n+1}\right)v + \int_{\Gamma_{T,N}} k_\mathrm{eff}^n N_T^{n+1} v + \int_{\Gamma_{\phi,N}} L_v^n\delta_\mathrm{p}^n\nabla p_\mathrm{sat}^n N_\phi^{n+1} v$

Doing the same process with a test function $q$ associated to $\phi$, the equation $\ref{eq:phi}$ becomes :

\[d(\phi^{n+1},q)=g(q)\]

with :

$d(\phi^{n+1},q) = \displaystyle\int_\Omega\xi^n\frac{\phi^{n+1}}{\Delta t}q + \int_\Omega\xi^nD_\mathrm{w}^n\nabla\phi^{n+1}\nabla q + \int_\Omega\delta_\mathrm{p}^n\nabla\phi^{n+1}p_\mathrm{sat}^n\nabla q + \int_\Omega\delta_\mathrm{p}^n\phi^{n+1}\nabla p_\mathrm{sat}^{n}\nabla q - \int_{\Gamma_{N,\phi}}\left(\delta_\mathrm{p}^n\phi^{n+1}\nabla p_\mathrm{sat}^n\cdot\mathbf{n}\right)q$
$g(q) = \displaystyle\int_\Omega\left(\xi^n\dfrac{\phi^n}{\Delta t} + G^{n+1}\right)q + \int_{\Gamma_{N,\phi}}\left(\xi^nD_\mathrm{w}^n+\delta_\mathrm{p}^np_\mathrm{sat}^n\right)N_\phi^{n+1}q$

2. Model properties

As explained in the Feel Tutorial Dev, the same set of equations allow solving many different problems. That is why we define a model in Feel.

Models are represented by a class ModelProperties which contains sub-classes corresponding to the different sections of a JSON file. For example, the following illustrates the definition of boundary conditions, for the process is quite similar for the other sections (initial conditions…)

This object is stored as a member of the class HeatMoisture, named props_.

3. Example of usage for boundary conditions

"BoundaryConditions":
{
    "field":                    (1)
    {
        "BC_type":              (2)
        {
            "marker":           (3)
            {
                "expr":"value"  (4)
            }
        }
    }

1	`field` is the field on what the boundary condition is applied (`phi` for relative humidity or `T` for heat)
2	`BC_type` is the type of boundary condition applied to the field (`Dirichlet` or `Neumann`)
3	`marker` is the name of the physical entity in the geo file, on which the condition is applied
4	`value` is the expression of the value of the condition applied to the field. For Dirichlet condition, it is equal to the value above it, and for Neumann condition, it is the value of $f$ such as $\frac{\partial \phi}{\partial n}=f$.

The boundary conditions are used in the code at the definition of the forms, in the functions update of the two classes (see section 5.3 of the project report)

We distinguich three type of physics used in the model :

heat, for heat transfer equations
moisture, for moisture transfer equations
heat-moisture, for the coupled heat and moisture transfer equations

We define the physics in the JSON file in a section Materials, on associated markers, allowing to combine physics in a more complex geometry :

"Materials":
{
    "Omega":
    {
        "physics":["moisture"],
        "markers":["Omega"],
    }
}

Let’s focus on the implementation in C++ of the boundary conditions. Here is the (compacted) code of the update function in the class Heat.

void update( double t, PhiT const& phi, double dt )
{
    this->a_.zero() ;
    this->l_.zero() ;

    for ( auto const& [name, mat] : this->props_.materials() ) (1)
    {
        if ( mat.hasPhysics( "heat" ) || mat.hasPhysics( "heat-moisture" ) ) (2)
        {
            // don't depend on boundary conditions
            this->a_ += integrate( _range = markedelements( this->mesh(), mat.meshMarkers() ),
                                   _expr = ... );
            this->l_ += integrate( _range = markedelements( this->mesh(), mat.meshMarkers() ),
                                   _expr = ... );

            if ( mat.hasPhysics( "moisture" ) || mat.hasPhysics( "heat-moisture" ) ) (2)
            {
                this->l_ += integrate( _range = markedelements( this->mesh(), mat.meshMarkers() ),
                                       _expr = ... );
            }

            // depend on the boundary conditions
            for (auto const& [bcid, bc] : this->props_.boundaryConditions2().flatten()) (3)
            {
                if (bcid.type() == "Neumann")
                {
                    if (bcid.field() == "T")    // for heat
                        this->l_ += integrate( _range = markedfaces( this->mesh(), bc.markers() ), _expr = ... );
                    if ( bcid.field() == "phi" && (mat.hasPhysics( "moisture" ) || mat.hasPhysics( "heat-moisture" )) )  // for moisture
                        this->l_ += integrate( _range = markedfaces( this->mesh(), bc.markers() ), _expr = .... );
                }
            }

            for (auto const& [bcid, bc] : this->props_.boundaryConditions2().flatten()) (4)
            {
                if (bcid.type() == "Dirichlet" && bcid.field() == "T")
                {
                    auto Ts = expr( bc.expression() ) ;
                    Ts.setParameterValues({"t",t}) ;
                    this->a_ += on( _range=markedfaces(this->mesh(), bc.markers()), _rhs=this->l_, _element=T_, _expr=Ts) ;
                }

            }
        }
    }
}

1	We begin by making a loop on all the materials.
2	If this material has the physic `heat` (or `heat-moisture`), we apply the conditions (and idem for the `moisture` term).
3	We make a loop on all the boundary conditions, adding in a first time the terms of Neumann conditions
4	Finally, we make a second loop on the boundary conditions, because we have to finish by adding the Dirichlet conditions.

4. Difficulties encountered

With all those modifications, the application works on the first try on the case drying-layer involving only the moisture transport process, but didn’t work on moisture-uptake : after a few steps, the linear solver of Feel++ didn’t converge, and therefore the results were totally erroneous !

When we try to compare the values of the coefficients obtained from this simulation and the ones obtained on the application when it was working (without the json file, it correspond to this this state in the history), we can see that problems appear atfer the first step. The main hypothesis to explain such an error is that the code struggles to handle some parameters because some expression depend on other expressions : for instance, K depends on w, and w depends on phi. A Feel++ fonctionnality allowing refactoring is under development and is currently available in the branch features/refactoring of the repository feelpp.

There were also some typos in the code such as missing parathensis. Once these problems solved, everything worked well !

5. Documentation of the application

The documentation of the application feelpp_hm_heat_moisture can be found here.

6. Test on several physics

Here, we study an imaginary case, with no real correspondance in physic. It is a mix between the two previous cases, moisture-uptake and drying-layer. The geometry is represented on figure 1 : on the left domain $\Omega_M$, the physic moisture of the drying-layer case is applied, and on the right one $\Omega_{HM}$, the physic heat-moisture is applied. Moreover, we a Dirichlet condition is added on the edge $\Gamma_\mathrm{mid}$ with $\phi=0.95$.

Figure 1. Geometry of the case

In figure 2 is represented the water content after a simulation of 7 days.

Figure 2. Water content (kg/m³)

Even though the results are not relevant in physics, we figure that they are coherent.

Application feelpp_hm_heat_moisture devlopment