An alternative equation-based model in COMSOL Multiphysics® for bentonite re-saturation

Michael Kröhn1, Lars Fromme2
1GRS gGmbH, Brunswick, Lower Saxony, Germany
2Bielefeld University of Applied Sciences, Faculty of Engineering and Mathematics, Bielefeld, Germany
Published in 2020

Bentonite is a versatile material that is, among other things, envisaged in designs for radioactive waste repositories to protect the waste canisters against groundwater. The thermo-hydraulic-mechanically coupled process of bentonite saturation is commonly based on two-phase flow and a mechanical stress formulation like the Barcelona Basic Model for unsaturated clayey soils. These features are also implemented in an example given in the COMSOL application library. As an alternative to the numerically demanding THM-formulations, a thermo-hydraulic saturation model for confined conditions (as expected in a repository) has been developed and experimentally realised with a one-dimensional formulation in the FORTRAN-code VIPER. In order to enhance the inherent limited range of possible applications of this code, the underlying partial differential equations are presently transferred to COMSOL Multiphysics® in form of an equation-based model. In the most simple, isothermal form, a partial differential equation has to be solved that looks like Fick’s second law for vapour diffusion which is non-linear, though, because the diffusion coefficient depends on the primary variable, the vapour partial density. Thus, only the mathematics interface was used, but couplings with other physic interfaces like the temperature interface should be implemented at a later stage. This balance equation was implemented using the coefficient form of the PDE interface. It was linearized by taking the solution of the previous timestep to calculate the diffusion coefficient for the actual timestep. As explained in the COMSOL blog, this required the previous solution operator. Additionally to this operator, also the Domain ODEs and DAEs interface were needed. In this interface, the “nojac”-command was used to write the saved value from the last solution into a new variable. At this stage the model input was set to match earlier calculations with code VIPER. These concern a 1D cross-section of a bentonite buffer that is wetted from one end and closed to flow at the other end (results: see Fig. 1). Validation of the new implementation in COMSOL was based on a comparison of the transient water content distributions calculated by both models. Based on this model, a first simple application was developed to provide a reference for subsequent modelling of multi-dimensional problems of the same type. A 2D-model simulating basically the 1D-problem (cp. Fig.2) was successfully tested that way. Setting up multi-dimensional problems in general was easier than expected, because only minor adjustments were required. The ability of true three-dimensional simulations are demonstrated by changing the Dirichlet boundary condition on the entire face of a 3D-slab to a point injection in the middle of the face, showing in the beginning a half-spherical expansion of the bentonite wetting (see Fig. 3). A next step will be the implementation of other physic interfaces to do real multiphysics calculations. Of particular interest is the coupling to the heat transfer interface and later also to a suitable mechanics interface. Presented here are the first results with COMSOL Multiphysics® demonstrating the gain of the transfer from the one-dimensional code VIPER in terms of multi-dimensional applications.