Finite element modelling of an energy-geomembrane underground pumped hydroelectric energy storage system

The increasing need for energy storage technology has led to a massive interest in novel energy storage methods. The energy geomembrane system is such a novel energy storage method. The concept of the system is briefly introduced, and a holistic numerical model of the system is presented. The model uses advanced finite-element techniques to model the energy storage system using fluid cavity elements. The developed geomembrane energy system is modelled with different constitutive models to represent the soil behaviour: a linear elastic model, a nonlinear Mohr-Coulomb model, and a hypoplastic constitutive model. The consequences of these different models on the results are studied. Hereby, the focus is the first inflation and deflation cycle of the system.


Introduction
The global effort to reduce and minimize anthropogenic greenhouse gas emissions requires an increase in the use of renewable energy sources. These energy resources are mainly produced by wind and solar power so far. A challenging aspect is the large fluctuations this sources of energy production can have with respect to supply and demand. Therefore energy storage systems are of high importance to ensure the green transition to minimize greenhouse gas emissions, which is extremely important to limit the global average temperature rise to 1.5 K to reach the Paris agreement [1].
However, energy storage systems that are adaptable in size (from residential-scale to communityscale) and can be used decentrally must be developed. Different scales and principle energy storage technologies are commonly used or are in development, e.g. shallow and deep geothermal energy storage systems [2,3,4], electrical energy storage systems, the transformation of energy (Power-to-Gas), and mechanical energy storage systems (compressed air reservoirs, pumped hydro reservoirs) [5]. Some of these systems are technically demanding. In a country like Denmark, the topography and geology are also essential factors, which limits the potential geotechnical energy storage solutions to a reasonable extent.
A system for storing volatile energy from renewable sources as gravitational potential energy is proposed by the company AquaNamic [6]. The novel system is called Underground Pumped Hydroelectric Energy Storage, abbreviated UPHS. One of the key advantages of this system is that it can be used in relatively flat topological areas found, for example, in Denmark, North Germany and the Netherlands.
The system operates by utilizing excess energy to store gravitational potential energy by inflating a subsurface reservoir with water from a nearby pond. The subsurface reservoir is lined with a geomembrane to form a bubble and the mass of the overburden soil is lifted as the geomembrane bubble is inflated. The inlets and outlets of fluid are controlled by a pumping system and turbine generators, similarly to the concept used in traditional pumped hydroelectric energy storage (PHES).
Within each cycle of inflation and deflation of the geomembrane bubble, mechanical energy dissipates due to different phenomena at the macro-scale (failure in the soil mass), at the micro-scale (intergranular friction), and within the interface between the soil and the geomembrane (see Figure 1).
Due to the novelty of this energy storage systems, this system is studied numerically to design and analyse the potentials for such a system. Therefore a holistic finite element model has been presented earlier in [7]. It was shown that the potential using fluid-cavity elements for holistically modelling the whole system is appropriate. In this model, only the hypoplastic constitutive model developed by v. Wolffersdorff [8] with the intergranular strain extension [9] has been used.
In this paper, the focus will be on the simulation of the energy system by using three different constitutive models for the soil behaviour, namely the elastic, the nonlinear Mohr-Coulomb model and the hypoplastic model. In addition, the different influences on the shear stress and the normal stress at the level of the geomembrane bubble are studied in greater detail. In the following section (2), the three constitutive models and the boundary value problem are introduced. Section 3.1 presents the results for the constitutive soil models and Section 3.2 the different results for the geomembranesoil interaction. Section 4 discussed the different results and concludes on the findings of this study.

Boundary value problem and constitutive models
The Energy membrane-UPHS (EM-UPHS) was modelled using the commercial finite element code ABAQUS (Dassault Systemes, version 2018). The numerical model was generated as an axisymmetric model considering the symmetry of the problem. ( Figure  2). The geometry of the boundary value problem is given in Table 1. The boundary problem itself has been generated using a parameterized Python script considering the design character of the energy storage problem. The modelling sequences are: Application of the geostatic stress field, followed by an initial first inflation and deflation cycle of the geomembrane which is modelled using the fluid cavity elements linking mechanism.
Besides the challenge of modeling the cyclic behaviour of the soil and ensuring the numerical convergence of the boundary value problem, the accurate modelling of the pumping mechanism was important. The pumping mechanism was modelled using fluid cavity elements for holistic modelling of the whole boundary problem. These fluid cavity elements were coupled with an auxiliary fluid container by a fluid link element, which can transfer fluid between the cavities. A schematic illustration is shown in Figure 3. More details about the numerical representation and the function of the pumping mechanism can be found in [7]. The underlying soil was modelled as stiff soil considering the compaction during the construction of the EM-UPHS. This hypothesis will be validated during a small scale field demonstrator project. The soilstructure interface between the geomembrane and the soil is modelled using a simple Mohr-Coulomb friction law. The frictional coefficient has been chosen to be quite small (μ = 0.15); this choice is justified since the energy dissipation is smaller than considering a larger friction coefficient at the soil-structure interface. Using such a small interface friction coefficient is a conservative choice because a higher coefficient will result in more stored energy. The numerical algorithm for the membrane-soil interaction was the mortar method in which at the soil surface, the slave is defined, and for the membrane, the master surface is applied. This modelling technique has been proven to be efficient for soil-structure interface modelling [10] compared to   more classical approaches such as, e.g., zero-thickness interface elements as [11,12]. The membrane itself is modelled using shell elements (SAX1) considering the thickness of the membrane. Finally, the soil was modelled using CAX4 elements with a reduced integration during some trial simulations; this choice of elements gave a satisfactory response in the simulation.
The constitutive models investigated for the soil were the following: the linear elastic (LE) model, the nonlinear Mohr-Coulomb (NLMC) model and the hypoplastic model of v. Wolffersdorff (HYPO) [8] with the intergranular strain extension [9]. The parameter sets can be found in Table 2. The overlying sand was modelled as fine sand in a dense state. The LE model uses a linear elastic relation. This model is an oversimplification for soils, but it can deliver some insights into the modelling processes. The model only requires two parameters: Young's modulus, E, and Poisson's ratio, ν.
Further improvement can be be achieved using the NLMC model uses of the elastoplastic framework and considers nonlinear Mohr-Coulomb envelopes [13]. Basically, this is a Mohr-Coulomb model where the friction angle is stress dependent, which is an important feature when an analysis includes a large stress level range. The model includes linear elastic response and it is implemented into the finite element method by use of implicit integration for the stress update. The model requires eight parameters.
The nonlinear HYPO model utilizes the framework of hypoplasticity considering a rate-independent constitutive equation [14]. Von Wolffersdorff [8] has formulated the version used in these simulations. This model is also using the intergranular strain concept after Niemunis & Herle [9]. More details about hypoplastic and modelling soils using hypoplasticity can be found, e.g., in [15]. In general, hypoplasticity has been used in many geotechnical engineering applications, e.g., [16,17]. The model requires 13 model parameters for the hypoplastic model (8 parameters) and the intergranular strain concept (5 parameters). In the following section, the results of the applications of the different constitutive models will be shown.

Impact of constitutive models on deformation pattern of the overburden soil
The modelling of the EM-UPHS was done with different constitutive models to prove the general concept and compare the capabilities for a further study of the energy storage system. Therefore, the analysis focused on the initial inflation and deflation cycle of the energy storage system. Figure 4 and Figure 5 show the responses of the EM-UPHS system with respect to the vertical and horizontal displacements during inflation and deflation respectively. Results are shown for the first inflation and deflation cycle, (ncy = 0.5 indicates full inflation and ncy = 1.0 a return to the inital deflated volumen of the geomembrane). It should be noticed that the membrane is never fully deflated, which is close to the planned operational stages of the energy storage system. The deformation pattern for the NLMC and the HYPO model are similar. In contrast, the LE model predicts smaller vertical and horizontal displacements than both of the other models.
From Figure 5 it can be observed that the simulations using the LE model are fully reversible, and it returns to the initial deformed configuration. On the   contrary, irreversible displacements in the soil for the inelastic models are observed. The hypoplastic model predicts slightly lower displacements in x-and ydirections than the NLMC model. This can be expected due to the ability to model the influence of the stress level (barotropy) and the influence of density (pycnotropy). In addition, the HYPO model includes the intergranular strain concept, which provides higher initial small strain stiffness and constitutive behaviour compared to the NLMC model. However, for the volumetric deformation of the upper soil part for each material model, the Mohr-Coulomb predicts a significant dilation during inflation and deflation. In contrast, the hypoplastic predicts a contraction during inflation and deflation, as seen in Figure 6. Notice that "calculation time" defines the different steps in the simulation of the BVP. For the elastic model, the volume more or less remains constant during inflation and deflation.

Impact of constitutive models on soil-membrane interface behaviour
In Figure 7, the contact pressures are shown. During inflation of the membrane, significant shear stresses around 30 kPa are observed close to the axis of symmetry. Furthermore, significant shear stresses of approximately same magnitude but the opposite sign is encountered at the flexion of the membrane, corresponding to upward displacement of the upper soil layer relative to the membrane. During deflation, a decrease in the interface shear stresses is observed.
The relative tangential slip of the three constitutive models is shown in Figure 8 for two points ηx = 0.08 (symmetry axis of the model) and ηx = 0.87 (membrane folding area). When gravity is applied, the linear elastic (LE) and Mohr-Coulomb (NLMC) model predicts approximately the same slip, which could indicate that initial yielding of the Mohr-Coulomb model has not occurred at this stage in the analysis, and the response, therefore, is elastic. The hypoplastic model predicts approximately two times larger slip than the Mohr-Coulomb model during application of gravity.
When inspecting the slip at a point at ηx = 0.08 after gravity is applied, the course of the slip curves for the two inelastic models is relatively identical, and the Mohr-Coulomb curve is merely a translation of the hypoplastic curve. That reflects the similarity between the two models, which was also observed for the deformation patterns in Figures 4 and 5. As was observed in Section 3.1, the deformation pattern of the elastic model differs significantly from the inelastic models, and this difference in the behaviour of the models is confirmed when inspecting the slip in the interfaces.
In Figure 9, the shear stress and normal stress in the interface are shown. The similarity in the behaviour of the Mohr-Coulomb and the hypoplastic model is not as apparent, as seen from the slip and overall deformation pattern. However, the Mohr-Coulomb response and the hypoplastic response are still closely coinciding. When inspecting the outermost point, more substantial shear stresses, i.e. more negative shear stresses, are predicted by the Mohr-Coulomb, which is also a consequence of the larger normal pressure on the interface predicted by the Mohr-Coulomb model. The elastic model is significantly different from the inelastic models. When compared with the Mohr-Coulomb model and the hypoplastic model, the elastic model underestimates the shear stress for the outermost point. Whereas the central point (ηx = 0.08) shows identical results for all three different models.

Conclusion & Discussion
The presented study concerned the finite element modelling of an EM-UPHS energy storage system with three different constitutive models. The pumping mechanism was realized using fluid cavity elements. By this procedure, the inflation and deflation could be modelled successfully [7]. This modelling strategy was used to examine the influence of the choice of the constitutive model on the simulated behaviour of the energy storage system. In the following paragraphs, the presented results are discussed in detail.
Linear elastic model (LE). The linear elastic model has the advantage of being simple and easy to calibrate. This simple constitutive model is very robust and less sensitive to the choice of parameters. The model provides a reasonable first estimate of the stress distribution in the soil. Due to the simplicity of the model, it is unable to capture the barotropic, pycnotropic and critical state behaviour, which soils possess. The stiffness of the soil and Poisson's ratio vary with confining pressure, and the model does not capture this effect. The model is fully reversible, i.e. no dissipation of energy and accumulation displacements and strains occur. This is far from the actual soil behaviour since strain energy dissipates during a hysteretic loading. The fact that the model has no tensile limit results in a poor prediction regarding the deformation behaviour. A poor prediction of the deformation behaviour leads to a poor prediction of the stored energy.
Elastoplastic Mohr-Coulomb model with nonlinear failure envelope (NLMC). This model has the advantage of having a failure envelope which is calibrated for peak states of the soil for different levels of confining pressure and, as a result of this, acknowledges that the friction angle of soil depends on the stress level. The model has a high convergence rate. This is probably due to the choice of parameters in this study, resulting in associated plastic flow, resulting in a quick and stable analysis. When inspecting the deformation behaviour, the model shows a consistent, realistic deformation behaviour, which is inelastic and therefore, the model also predicts plastic dissipation. When inspecting the dilative behaviour of the model, excessive unrealistic dilation is predicted. However, calibration of a plastic potential, which is determined from deformation parameters, could probably improve the dilative behaviour. If the cohesion in the model was reduced to a value below 3 kPa, the models proved to lack computational robustness. Furthermore, the model is sensitive to changes in the parameters used to define the failure envelope.
Hypoplastic model with intergranular strain (HYPO). This model is able to capture pycnotropy and barotropy most successfully of the models examined, as it utilizes state variables to determine the soil response. It can simulate reasonable dilative and contractile behaviour and predicts realistic deformation behaviour. The parameters which should be defined for the model are invariant of stress states and the denseness of the soil and characteristic of the material. The model proved to be sufficiently robust for a small value of the cohesion (1 kPa). The model shows a highly nonlinear irreversible behaviour, which allows for simulating strain energy dissipation in the soil. A shortcoming of the model is the rather complex formulation. Due to the high nonlinearity of the model, the convergence rate is slow. However, the model still converged. Some of the model parameters were quite sensitive, and often require slight adaptations for numerical convergence, which seems to be reasonable for a quite complex nonlinear boundary value problems, with a lot of different contact considerations, as the one presented in this study.
The HYPO and NLMC models have shown that a realistic constitutive model is needed to model all the different phenomena that potentially will influence the energy storage system. An aspect which has not been studied in detail so far is the more advanced modelling of the geomembrane-soil interface. This will be done in future studies in connection with a comparison with small-field scale tests to verify the whole system in a small demonstrator scale.