Maximum entropy based finite element analysis of porous media

Abstract

The maximum entropy theory has been used in a wide variety of physical, mathematical and engineering applications in the past few years. However, its application in numerical methods, especially in developing new shape functions, has attracted much interest in recent years. These shape functions possess the potential for performing better than the conventional basis functions in problems with randomly generated coarse meshes. In this paper, the maximum entropy theory is adopted to spatially discretize the deformation variable of the governing coupled equations of porous media. This is in line with the well-known fact that higher-order shape functions can provide more stable solutions in porous problems. Some of the benchmark problems in deformable porous media are solved with the developed approach and the results are compared with available references.

Keywords

maximum entropy FEM / fully coupled multi-phase system / porous media

Introduction

Studying porous materials to characterize their mechanical behavior in different conditions is essential in many engineering applications, including groundwater flow, wells, oil reservoirs, biological tissues, nano-scale materials, etc.

Movement of fluid in a porous medium can cause deformation in its solid structure. While this phenomenon may be less important in many engineering problems, deformation of solid structures is crucially important in certain civil engineering applications to ensure that solid foundations remain stable, as consolidation of foundation and settlement of structures can lead to severe damages of buildings.

Analytical solution of the multiphase flow in porous media may be applicable in very simple problems using some simplified assumptions to reduce the complexity of non-linear equations. In general, however the governing coupled differential equations for deformable porous systems have been solved numerically using, for instance, the finite difference method [¹–³], the finite volume method [⁴–⁶] and discontinuous Galerkin methods [⁷,⁸], mostly used for modelling of flow and transport of miscible and non-miscible fluids in porous media. For coupled deformation analysis, the finite element method [⁹–¹⁵], meshless techniques [¹⁶–¹⁹] and XFEM [²⁰–²²] have been frequently used for continuum applications, discontinuities and fracture problems. Moreover, application of multiscale homogenization modelling in porous media has been widely increased in recent years [²³–²⁶].

In a pioneering work, Beltzer evaluated the complexity of finite element method with the maximum entropy [²⁷]. Sukumar applied the Shannon’s entropy [²⁸] to obtain the minimum–biased interpolants on polygonal domains to construct polygonal shape functions [²⁹] and Arroyo and Ortiz modified the local maximum entropy basis for numerical solution of PDEs [³⁰]. Millán et al. used the Galerkin method with cell-based maximum entropy basis functions for PDEs to construct a smooth response with a good control on unstructured meshes [³¹].

Ortiz et al. showed the robustness of the maximum entropy meshless method with a tetrahedral background mesh for integration of formulations of incompressible problems and locking-free small strain elasticity [³²,³³]. Quaranta et al. used the maximum entropy basis function to analyze nonlinear reinforced concrete shear walls [³⁴]. Ullah et al. proposed an adaptive finite element-maximum entropy method, to refine the initial finite element model with the maximum entropy approach [³⁵]. Moreover, the method was adopted for applied fracture problems [³⁶–³⁸].Wu et al. applied the adaptive method for material and geometrical nonlinearities to solve the convection-diffusion problem [³⁹].

Application of the maximum entropy meshless method in incremental small strain elastoplastic analysis of geotechnical models was reported by Kardani et al. [⁴⁰] and Nazem et al. employed this method for consolidation of porous media and investigated its robustness and stability [⁴¹]. Navas et al. proposed a meshfree porous media model based on the maximum entropy principle and validated a set of benchmark steady seepage problems as well as static and dynamic consolidation with the B-bar algorithm to prevent the locking of the fluid phase [⁴²–⁴⁴]. Also, Zakrzewski et al. used this method to simulate an undrained layer of soil with large deformations [⁴⁵].

In numerical modelling of thermo-hydro-mechanical coupled problems, which the major variables are temperature, pressure and displacement, it is generally accepted to adopt a higher order interpolation scheme for displacement [⁴⁶]. In the conventional finite element method, this is usually performed by higher order Lagrangian basis shape functions, which require more nodes per element. The maximum entropy shape functions, however, can provide higher order continuity with the same number of nodes. In this paper, the work of Sukumar [²⁹] is extended to porous media and its performance is assessed with available reference analytical and numerical results.

After this introduction, the maximum entropy shape function is presented in Section ‎2. Then, the governing equations of porous media are described in Section ‎3. Section 4 is dedicated to solution and discussion of the benchmark problems of porous media. The concluding remarks are presented in Section 5.

Maximum entropy shape functions

The concept of entropy in the information theory was defined by Shannon [²⁸] for measuring the uncertainty of data or the insufficiency of knowledge. The maximum entropy concept was then introduced by Jaynes [⁴⁷] based on the Shannon theory [²⁸] as the least biased statistical inference of an event occurrence. He discussed the relation between the maximum entropy and other spectral methods and concluded that this method could have an optimal result where prior information was available. The concept of maximum entropy has been applied to a wide range of applications, from atomic and molecular problems and nuclear physics [²⁸] to image processing [⁴⁸] and economics [⁴⁹], as a good mean for fields with insufficient data.

Shannon proposed a robust method for making decision and predicting occurrence of an event in the case of uncertainty in data interpolation problems [²⁸].

(1)H(p_{\mathrm{1}},p_{\mathrm{2}},\mathrm{...},p_n)=-k_B{\displaystyle \sum _{i=\mathrm{1}}^n{p}_i\log p_i},

where p_i represents the probability of phase case i to happen and k_B is the Boltzmann constant. A higher entropy indicates the most probable state, which is associated with more disorder [⁵⁰].

H(p_{\mathrm{1}},p_{\mathrm{2}},\mathrm{...},p_n)should be a continuous monotonically increasing function for the case of equal probabilities. This condition can well be satisfied by the logarithmic form. The entropy of a continuous distribution, as defined in the information technology systems, can be written as Ref. [²⁸]:

(2)H=-{\displaystyle {\int }_{-\infty }^{+\infty }p(x)}\log p(x)dx\mathrm{.}

H is a good measure of uncertainty and Jaynes showed that maximizing the defined entropy would lead to the most unbiased state [⁵¹]. The term p(x) is the density distribution function and states the probability of occurrence. It is noted that k_B is considered unity in the informational entropy and is the key point of difference with the thermodynamic entropy defined in Eq. (1). Similarly, the entropy for a discrete set of probabilities of an event i is written in the form of Eq. (3)

(3)H=-{\displaystyle \sum _{i=\mathrm{1}}^n{p}_i\log p_i}\mathrm{.}

It is clear that the summation on non-negative probabilities is always one; {\displaystyle \sum _{i=\mathrm{1}}^n{p}_i}=\mathrm{1}, which represents the concept of partition of unity. This fact can be considered as a constraint to the maximizing problem. Also, there may be extra information, which should be considered in the form of following constraints in terms of x and y for each state.

(4){\displaystyle \sum _{i=\mathrm{1}}^n{p}_i}=\mathrm{1},

(5)x={\displaystyle \sum _{i=\mathrm{1}}^n{p}_i{x}_i},

(6)y={\displaystyle \sum _{i=\mathrm{1}}^n{p}_i{y}_i},

where x and y are the coordinates of a point inside the element andx_iandy_i are the nodal coordinates. The Lagrange multipliers method is employed to impose the mentioned constraints. Therefore, the functional L is defined as:

(7)L=-{\displaystyle \sum }_{i=\mathrm{1}}^n{p}_i\log p_i+(\alpha -\mathrm{1})\left(\mathrm{1}-{\displaystyle \sum }_{i=\mathrm{1}}^n{p}_i\right)+\beta \left(x-{\displaystyle \sum }_{i=\mathrm{1}}^n{p}_i{x}_i\right)+\gamma \left(y-{\displaystyle \sum }_{i=\mathrm{1}}^n{p}_i{y}_i\right)\mathrm{.}

Differentiating the functional L with respect to p_i, maximizes the entropy and fits the constraints. The Lagrange multiplier\alpha -1 is assumed in this form to simplify the manipulation process:

(8){\displaystyle \frac{\partial L}{\partial p_i}}=-\mathrm{1}-\log p_i-(\alpha -\mathrm{1})-\beta (x_i)-\gamma (y_i)=\mathrm{0},\mathrm{f}\mathrm{o}\mathrm{r}i=1,2,\mathrm{3...}

Rewriting Eq. (8) leads to:

(9)\log p_i=-\alpha -\beta x_i-\gamma y_i\to p_i=e^{-\alpha -\beta x_i-\gamma y_i}

The calculatedp_i corresponds to the maximum uncertainty, defined in Eq. (3), along with constraints of Eqs. (4) to (6). The partition of unity Eq. (4) is used to obtain:

(10){\displaystyle \sum }_{i=\mathrm{1}}^n{p}_i=e^{-\alpha }{\displaystyle \sum }_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}=\mathrm{1}\to e^{-\alpha }={\displaystyle \frac{\mathrm{1}}{{{\displaystyle \sum }}_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}}}\mathrm{.}

Substituting the last term of Eq. (10) into Eq. (9), the weighed form of probability is obtained,

(11)p_i={\displaystyle \frac{e^{-\beta x_i-\gamma y_i}}{{{\displaystyle \sum }}_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}}}\mathrm{.}

Here, the two remaining coefficients are calculated by multiplying p_i to e^{\alpha }{x}_i and e^{\alpha }{y}_i, respectively, and a summation on the probabilities leads to:

(12)e^{\alpha }{x}_i{p}_i=e^{-\beta x_i-\gamma y_i}{x}_i\Rightarrow e^{\alpha }{\displaystyle \sum }_{i=\mathrm{1}}^n{x}_i{p}_i={\displaystyle \sum }_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}{x}_i\stackrel{x={\displaystyle \sum }_{i=\mathrm{1}}^n{x}_i{p}_i}{\Rightarrow }{\displaystyle \frac{{{\displaystyle \sum }}_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}{x}_i}{e^{\alpha }}}-x=\mathrm{0},

(13)e^{\alpha }{y}_i{p}_i=e^{-\beta x_i-\gamma y_i}{y}_i\Rightarrow e^{\alpha }{\displaystyle \sum }_{i=\mathrm{1}}^n{y}_i{p}_i={\displaystyle \sum }_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}{y}_i\stackrel{y={\displaystyle \sum }_{i=\mathrm{1}}^n{y}_i{p}_i}{\Rightarrow }{\displaystyle \frac{{{\displaystyle \sum }}_{i=\mathrm{1}}^n{e}^{-\beta x_i-\gamma y_i}{y}_i}{e^{\alpha }}}-y=\mathrm{0.}

b and g are obtained by solving the two Eqs. (12) and (13) simultaneously. This step becomes highly complex if more constraints are involved. The numerical method developed by Sukumar [²⁹] is adopted for solving these set of equations. The calculated probabilities can be considered as finite element shape functions, because they possess the necessary condition of partition of unity.

The exponential form of obtained shape functions leads to higher order of continuity. The developed shape functions {{\displaystyle \mathit{N}}}_u^H (in terms of {{\displaystyle p}}_i) also can be used for interpolating the displacement field:

(14)u={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{N}}_u^H{\widehat{u}}_i},{\mathit{N}}_u^H=p_i\mathrm{.}

It is noted that the entropy based shape functions become bi-linear for a rectangular finite element, whereas they become highly non-linear in terms of exponential functions for general distorted elements [²⁹]. Figure 1 shows the generated shape functions on a typical distorted element, with coordinates of (0,0), (1,0.5), (0.5,1) and (1,1), showing a major difference in comparison with the conventional finite element shape functions.

Governing equations

Linear momentum, mass conservation (continuity equation) and energy conservation (enthalpy equation) govern the behavior of porous media. All the assumptions, the way the equations are derived and the finite element weak forms are well covered in many references [¹⁸,⁴⁶,⁵²,⁵³].

Beginning with the list of governing equations for an unsaturated multi-phases porous medium:

Linear momentum conservation:

(15){{\mathit{L}}_{\mathit{u}}}^{\mathrm{T}}\mathrm{.}\sigma +\rho \mathit{b}=0,

Linear elastic isotropic constitutive relation:

(16){\sigma }^{\prime \prime }=D_eϵ,

Darcy low:

(17)v^{\alpha s}={\displaystyle \frac{\mathit{k}{k}_{r\alpha }}{{\mu }_{\alpha }}}(-\nabla P_{\alpha }+{\rho }_{\alpha }\mathit{b}),

Mass conservation for phase α (no mass transfer):

(18){\displaystyle \frac{D}{Dt}}(n_{\alpha }{\rho }_{\alpha })+n_{\alpha }{\rho }_{\alpha }\nabla v_{\alpha }=\mathrm{0},

where vector bis the body force, \rhois the average density of system, \sigma, {\sigma }^{″}and \epsilonare the total stress, effective stress and strain tensors respectively, D_eis the elastic stiffness matrix, v^{\alpha s} is the relative velocity of phase a to solid, k is the intrinsic permeability, k_{r\alpha }is the relative permeability of phase a and {\mu }_{\alpha } is the viscosity of phase a.

The initial conditions are:

(19)u_{(t=0)}=0P_{w(t=0)}=P_{w\left(0\right)}{P}_{g(t=0)}=P_{g\left(0\right)},

The essential boundary conditions are:

(20)u=\overline{u}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_u,

(21)P_w=\overline{P_w}\mathrm{on}{\Gamma }_{P_w},

(22)P_g=\overline{P_g}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{P_g},

and the natural boundary conditions can be written as:

(23)\overline{t}=\sigma .n_{{\Gamma }_t}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_t,

(24){\overline{q}}_w=v_w.n_{{\Gamma }_w}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{q_w},

(25){\overline{q}}_g=v_g.n_{{\Gamma }_g}on{\mathrm{\Gamma }}_{q_g},

where \bar{u}, \overline{P_w}and \overline{P_g}are the predefined values of the independent variables on the essential boundary conditions {\Gamma }_u, {\Gamma }_{P_w}and {\Gamma }_{P_g}, respectively. Also,\bar{t}, {\bar{q}}_w and {\bar{q}}_gare the rate of external force and flux vectors on the natural boundary conditions {\Gamma }_t,{\Gamma }_{q_w}and {\Gamma }_{q_g}, respectively.

For an unsaturated porous medium, the following discretized set of coupled simultaneous equations are obtained [⁵⁴]:

(26)\left[\begin{array}{lll}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ C_{ws}& P_{ww}& C_{wg}\\ C_{gs}& C_{gw}& P_{gg}\end{array}\right]{\displaystyle \frac{d}{dt}}\left[\begin{array}{l}u\\ P_w\\ P_g\end{array}\right]+\left[\begin{array}{lll}{K}_e& -C_{sw}& -C_{sg}\\ \mathrm{0}& H_{ww}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& H_{gg}\end{array}\right]\left[\begin{array}{l}u\\ P_w\\ P_g\end{array}\right]=\left[\begin{array}{l}{F}_u\\ F_w\\ F_g\end{array}\right],

where the maximum entropy shape function N_u^His used for the displacement discretization, while the classical finite element shape function N_P is used for discretizing the pressure and:

(27)K_e={\displaystyle {\int }_{\Omega }{(B_u^H)}^T{D}_e(B_u^H)\mathrm{d}\mathrm{\Omega }},

(28)C_{sw}={\displaystyle {\int }_{\Omega }{(B_u^H)}^T\alpha m{S}_w{N}_P\mathrm{d}\mathrm{\Omega }},

(29)C_{sg}={\displaystyle {\int }_{\Omega }{(B_u^H)}^T\alpha m{S}_g{N}_P\mathrm{d}\mathrm{\Omega }},

(30)F_u={\displaystyle {\int }_{\Omega }{(N_u^H)}^{\mathrm{T}}\rho g\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\Gamma }{(N_u^H)}^{\mathrm{T}}\overline{t}\mathrm{d}\mathrm{\Gamma }}

with

(31)H_{ww}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{B}}_P}^{\mathrm{T}}}\frac{{\mathit{k}k}_{rw}}{{\mu }_w}{\mathit{B}}_P\mathrm{d}\mathrm{\Omega },

(32)C_{ws}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{N}}_P}^{\mathrm{T}}\alpha S_W{\mathit{m}}^{\mathrm{T}}\left({{\mathit{B}}_u}^H\right)}\mathrm{d}\mathrm{\Omega },

(33)P_{ww}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{N}}_P}^{\mathrm{T}}\left[\frac{S_w(\alpha -n)}{K_S}(S_w+P_c\frac{\partial S_w}{\partial P_c})+n\frac{S_w}{K_w}-n\frac{\partial S_w}{\partial P_c}\right]{\mathit{N}}_P}\mathrm{d}\mathrm{\Omega },

(34)C_{wg}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{N}}_P}^{\mathrm{T}}\left[\frac{S_w(\alpha -n)}{K_S}(S_g-P_c\frac{\partial S_w}{\partial P_c})+n\frac{\partial S_w}{\partial P_c}\right]{\mathit{N}}_P}\mathrm{d}\mathrm{\Omega },

(35)F_w={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{B}}_P}^{\mathrm{T}}}\frac{{\mathit{k}k}_{rw}}{{\mu }_w}{\rho }_{w}g\mathrm{d}\mathrm{\Omega }-{\displaystyle {\int }_{\mathrm{\Gamma }}{{\mathit{N}}_P}^{\mathrm{T}}}{\overline{q}}_w\mathrm{d}\mathrm{\Gamma },

and

(36)H_{gg}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{B}}_P}^{\mathrm{T}}}\frac{{\mathit{k}k}_{rg}}{{\mu }_g}{\mathit{B}}_P\mathrm{d}\mathrm{\Omega },

(37)C_{gs}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{N}}_P}^{\mathrm{T}}\alpha S_g{\mathit{m}}^{\mathrm{T}}\left({{\mathit{B}}_u}^H\right)}\mathrm{d}\mathrm{\Omega },

(38)C_{gw}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{N}}_P}^{\mathrm{T}}\left[\frac{S_g(\alpha -n)}{K_S}(S_w+P_c\frac{\partial S_w}{\partial P_c})+n\frac{\partial S_w}{\partial P_c}\right]{\mathit{N}}_P}\mathrm{d}\mathrm{\Omega },

(39)P_{gg}={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{N}}_P}^{\mathrm{T}}\left[\frac{S_g(\alpha -n)}{K_S}(S_g-P_c\frac{\partial S_w}{\partial P_c})-n\frac{\partial S_w}{\partial P_c}+n\frac{S_g}{K_g}\right]{\mathit{N}}_P}\mathrm{d}\mathrm{\Omega },

(40)F_g={\displaystyle {\int }_{\mathrm{\Omega }}{{\mathit{B}}_P}^{\mathrm{T}}}\frac{{\mathit{k}k}_{rg}}{{\mu }_g}{\rho }_{g}g\mathrm{d}\mathrm{\Omega }-{\displaystyle {\int }_{\mathrm{\Gamma }}{{\mathit{N}}_P}^{\mathrm{T}}}{\overline{q}}_g\mathrm{d}\mathrm{\Gamma },

where n is the porosity of material,S_w and S_g are the water and gas saturation, ratio a is Biot constant, \mathit{m}={\left[\begin{array}{lll}\mathrm{1}& \mathrm{1}& \mathrm{0}\end{array}\right]}^T, K_s,K_wand K_g are the bulk modulus of solid, water and gas respectively.

{\mathit{B}}_u^H={\mathit{L}}_u{\mathit{N}}_u^Hand {\mathit{B}}_P={\mathit{L}}_P{\mathit{N}}_P, where {\mathit{L}}_uand{\mathit{L}}_Pare the differential operators:

(41){\mathit{L}}_u=\left[\begin{array}{ll}{\displaystyle \frac{\partial }{\partial x}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\partial }{\partial y}}\\ {\displaystyle \frac{\partial }{\partial y}}& {\displaystyle \frac{\partial }{\partial x}}\end{array}\right],{\mathit{L}}_P=\left[\begin{array}{l}{\displaystyle \frac{\partial }{\partial x}}\\ {\displaystyle \frac{\partial }{\partial y}}\end{array}\right]\mathrm{.}

P_c is the capillary pressure and both P_c and {\displaystyle \frac{\partial S_w}{\partial P_c}} are chosen from the experimental relations [¹⁰,⁵⁵].

Eq. (28) can be rewritten as:

(42)B{\displaystyle \frac{dX}{dt}}+CX=F,

which is a nonlinear differential equation in time. This time-dependent set of equations is solved in time by the fully implicit method [⁴⁶]:

(43){\left[{\displaystyle \frac{B}{\Delta t}}+C\right]}_{n+\mathrm{1}}{X}_{n+\mathrm{1}}-{\left[{\displaystyle \frac{B}{\Delta t}}\right]}_{n+\mathrm{1}}{X}_n=F_{n+\mathrm{1}}\mathrm{.}

A Newton-Raphson scheme is adopted to solve the nonlinear set of Eq. (43) [⁴⁶].

Numerical simulations

In this section, benchmark porous media problems are solved with the maximum entropy shape functions in order to verify the developed algorithm and the results are compared with the available reference results.

Saturated porous medium

Two-dimensional saturated foundation under a strip loading

The first example is the consolidation of a layer of soil under a strip loading. A semi analytical solution of this example was derived by Brooker and Small [⁵⁶,⁵⁷] using the finite difference method and the Fourier transform of governing equations of saturated porous media.

The length of foundation is 12m and only the half of the problem is modeled, due to the symmetry of problem (Fig. 2). The layer is 1 m height and 6 m length and 1000 Pa magnitude strip loading is applied to the foundation between 0≤a≤1 m.This problem is now solved for two types of boundary conditions, defined in Table 1. The material properties are presented in Table 2.

The problem is solved for the total time t_{\mathrm{total}}=10,000 s and with the time step \Delta t=0.1 s. The results are compared with the analytical solution of Gibson et al. [⁵⁸] and numerical solution of Samimi and Pak [¹⁶]. Figure 3 and 4 present changes of surface settlement and water pressure, respectively, along the middle height of symmetric line in time. Clearly, they are in good agreement with the references results, even though a randomly generated mesh has been adopted.

Thermo-elastic porous medium

The second problem is about a saturated elastic porous medium under the heat loading. In a thermo saturated porous medium, all phases are locally in a state of thermodynamic equilibrium [⁴⁶]. With this assumption, the general form of enthalpy equation (without considering heat transfer) can be written as:

(44){(\rho C)}_{\mathrm{eff}}{\displaystyle \frac{\partial T}{\partial t}}+({\rho }_w{C}_w{v}_w+{\rho }_g{C}_g{v}_g).\nabla T-{\lambda }_{\mathrm{eff}}.{\nabla }^{\mathrm{2}}(T)=\mathrm{0},

with the effective heat capacity {(\rho C)}_{\mathrm{eff}}:

(45){(\rho C)}_{\mathrm{eff}}=(\mathrm{1}-n){\rho }_s{C}_s+n{\rho }_w{C}_w,

Also, the effective heat conductivity is:

(46){\lambda }_{\mathrm{eff}}=(\mathrm{1}-n){\lambda }_s+n{\lambda }_w,

where C_{\alpha } is the heat capacity and {\lambda }_{\alpha }is the heat conductivity of phase a.

This example, previously solved by Aboustit et al. [⁵⁹] and Lewis and Schrefler [⁴⁶], is comprised of an elastic soil column under a vertical load on the top of the column. A temperature increase of \Delta T=\mathrm{50} K is also applied on the surface of soil. Figure 5 shows the finite element mesh and part of boundary conditions. Four-node elements with the size of 25 mm×25 mm are used and each element has four Gauss points.

The boundary conditions are defined in Table 3, and the material properties are shown in Table 4. The time step is set to \Delta t=0.1 s and the total time is t_{\mathrm{total}}=10,000 s.

Figures 6, 7, and 8 depict the changes of water pressure, temperature and vertical displacement, respectively, in time at distances of 0.2 m, 1 m, and 3 m from the surface. The results are in agreement with the reference results. Figure 9 outlines the contours of variation of temperature, which shows how the temperature changes across the height of the specimen in time and towards a uniform temperature.

Unsaturated porous medium

In this section, an unsaturated porous medium test, known as the Liakopolous problem, is studied by two methods. The test performed on Del Monte sand by Liakopoulos [⁶⁰], is recognized as a benchmark problem in non-saturated porous media, as the sample becomes non-saturated due to water discharge from the bottom of sample. This experiment has been examined by several researchers numerically [¹¹,¹⁹,⁵²–⁵⁴,⁶¹–⁶⁷]. The experimental conditions can be assumed as a uniaxial sample column with 1 m height, which is filled by sand and then the water is poured from the top of sample, until the porous medium becomes saturated. Then, the inflow rate of water is ceased and the water is discharged from the bottom of sample, changing the state of porous medium to an unsaturated condition. The test is solved without and in the presence of air flow in the following sections. The material parameters are listed in Table 5.

A Porous medium with deactivated air flow

In this case, the independent unknowns are displacement and water pressure and the air pressure is assumed to remain constant at atmospheric pressure all over the sample. Hence, the capillary pressure isP_c=-P_w. In addition, it is assumed that after the inflow of water is ceased, the pressure at the top of sample remains constant at the atmospheric pressure. The finite element mesh of problem is shown in Figure 10.

The initial capillary pressure is zero at t=0 s in every nodes. Water pressure and displacement boundary conditions are defined in Table 6.

As illustrated in Figure 10, the problem is solved with a 25 cm×25 cmmesh and each element has 4 Gauss points. In order to calculate the saturation and water relative permeability equations in every time steps and every Gauss points, the following empirical relation is adopted [⁵⁴]:

S_w=\mathrm{1}-\mathrm{0.10152}{\left({\displaystyle \frac{P_c}{{\rho }_{w}g}}\right)}^{\mathrm{2.4279}}\mathrm{for}{S}_w\ge \mathrm{0.91},

k_{rw}=\mathrm{1}-\mathrm{2.207}{(\mathrm{1}-S_w)}^{\mathrm{1.0121}}\mathrm{.}

The problem is solved by Δt=20 s for the total time of two hours and the results of water pressure, capillary pressure, saturation and vertical displacement are presented in Figs. 11‒14, respectively. The results illustrate variations of each variable in different times versus the height of sample. Clearly, the results are perfectly compatible with the reference results.

A porous medium with activated air flow

The problem of 4.2.1 is solved again with the assumption of change of air pressure along the column. Therefore, the air pressure is assumed as an independent variable in the finite element formulation. The air pressure at the top and bottom of sample is assumed constant (atmospheric pressure). Hence, the capillary pressure is equal to the difference of air and water pressures: P_c=P_{\mathrm{air}}-P_w. The material properties are defined in Table 5. Also, the adopted finite element mesh is depicted in Fig. 15.

The capillary pressure is zero at each node due to imposition of initial conditions. Displacement, water pressure and air pressure boundary conditions of problem 4.2.2 are defined in Table 7.

The elements are similar to the previous example and the time step and total time are \Delta t=0.1 s andt_{\mathrm{total}}=7200 s, respectively. For calculation of water pressure and relative permeability of water, Eqs. (47) and (48) are used. For the air relative permeability, the Brooks and Corey relation is utilized [⁵⁵]:

(49)s_e=(s_w-s_{wr})/(\mathrm{1}-s_{wr}),

(50)s_e=(s_w-s_{wr})/(\mathrm{1}-s_{wr}),

wheres_e, s_{wr}=0.2 and \lambda=3 are the effective water saturation, the residual water saturation and the pore size distribution index, respectively. The results for 5, 10, 20, 30, 60, and 120 minutes are depicted in Figs. 16‒20.

Finally, the outflow rate of both methods are obtained in time and compared with the laboratory data in Fig. 21 to examine the best performance.

Conclusions

In this paper, the recently developed basis of maximum entropy concept is used to construct the shape functions of finite element method, to study the behavior of porous problems. The shape functions have the capability of solving the problem even in distorted meshes and are capable of obtaining higher orders of continuity with the same degrees of freedom.

Some of the major problems in saturated, thermo-saturated and unsaturated porous media have been studied with this method. The results are in good agreement with available references.

The method has been used in the field of porous media and coupled problems. Nevertheless, it is potentially attractive in other complex problems such as discontinuities and fracture. On the other hand, the complexity of solving equations and high expenses of constructing the shape functions, should be realistically examined for general engineering applications.

References

[1] Touma J, Vauclin M. Experimental and numerical analysis of two-phase infiltration in a partially saturated soil. Transport in Porous Media, 1986, 1(1): 27–55

[2] Faust C R, Guswa J H, Mercer J W. Simulation of three-dimensional flow of immiscible fluids within and below the unsaturated zone. Water Resources Research, 1989, 25(12): 2449–2464

[3] Ataie-Ashtiani B, Raeesi-Ardekani D. Comparison of numerical formulations for two-phase flow in porous media. Geotechnical and Geological Engineering, 2010, 28(4): 373–389

[4] Durlofsky L J. A triangle based mixed finite element—finite volume technique for modeling two phase flow through porous media. Journal of Computational Physics, 1993, 105(2): 252–266

[5] Forsyth P A, Wu Y, Pruess K. Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media. Advances in Water Resources, 1995, 18(1): 25–38

[6] Jenny P, Lee S H, Tchelepi H A. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. Multiscale Modeling & Simulation, 2005, 3(1): 50–64

[7] Klieber W, Rivière B. Adaptive simulations of two-phase flow by discontinuous Galerkin methods. Computer Methods in Applied Mechanics and Engineering, 2006, 196(1‒3): 404–419

[8] Epshteyn Y, Rivière B. Fully implicit discontinuous finite element methods for two-phase flow. Applied Numerical Mathematics, 2007, 57(4): 383–401

[9] .Li X O, Zienkiewicz . Multiphase flow in deforming porous media and finite element solutions. Computers & Structures, 1992, 45(2): 211–227

[10] Rahman N A, Lewis R W. Finite element modelling of multiphase immiscible flow in deforming porous media for subsurface systems. Computers and Geotechnics, 1999, 24(1): 41–63

[11] Laloui L, Klubertanz G, Vulliet L. Solid–liquid–air coupling in multiphase porous media. International Journal for Numerical and Analytical Methods in Geomechanics, 2003, 27(3): 183–206

[12] Oettl G, Stark R, Hofstetter G. Numerical simulation of geotechnical problems based on a multi-phase finite element approach. Computers and Geotechnics, 2004, 31(8): 643–664

[13] Stelzer R, Hofstetter G. Adaptive finite element analysis of multi-phase problems in geotechnics. Computers and Geotechnics, 2005, 32(6): 458–481

[14] Callari C, Abati A. Finite element methods for unsaturated porous solids and their application to dam engineering problems. Computers & Structures, 2009, 87(7‒8): 485–501

[15] Nguyen V P, Lian H, Rabczuk T, Bordas S. Modelling hydraulic fractures in porous media using flow cohesive interface elements. Engineering Geology, 2017, 225: 68–82

[16] Samimi S, Pak A. Three-dimensional simulation of fully coupled hydro-mechanical behavior of saturated porous media using Element Free Galerkin (EFG) method. Computers and Geotechnics, 2012, 46: 75–83

[17] Goudarzi M, Mohammadi S. Weak discontinuity in porous media: an enriched EFG method for fully coupled layered porous media. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38(17): 1792–1822

[18] Goudarzi M, Mohammadi S. Analysis of cohesive cracking in saturated porous media using an extrinsically enriched EFG method. Computers and Geotechnics, 2015, 63: 183–198

[19] Samimi S, Pak A. A three-dimensional mesh-free model for analyzing multi-phase flow in deforming porous media. Meccanica, 2016, 51(3): 517–536

[20] Mohammadnejad T, Khoei A. Hydro-mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37(10): 1247–1279

[21] Goodarzi M, Mohammadi S, Jafari A. Numerical analysis of rock fracturing by gas pressure using the extended finite element method. Petroleum Science, 2015, 12(2): 304–315

[22] Mohammadnejad T, Khoei A R. An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elements in Analysis and Design, 2013, 73: 77–95

[23] Zhuang X, Wang Q, Zhu H. A 3D computational homogenization model for porous material and parameters identification. Computational Materials Science, 2015, 96: 536–548

[24] Zhu H, Wang Q, Zhuang X. A nonlinear semi-concurrent multiscale method for fractures. International Journal of Impact Engineering, 2016, 87: 65–82

[25] Bayesteh H, Mohammadi S. Micro-based enriched multiscale homogenization method for analysis of heterogeneous materials. International Journal of Solids and Structures, 2017, 125: 22–42

[26] Fatemi Dehaghani P, Hatefi Ardakani S, Bayesteh H, Mohammadi S. 3D hierarchical multiscale analysis of heterogeneous SMA based materials. International Journal of Solids and Structures, 2017, 118–119: 24–40

[27] Beltzer A I. Entropy characterization of finite elements. International Journal of Solids and Structures, 1996, 33(24): 3549–3560

[28] Shannon C E. Communication theory of secrecy systems. Bell Labs Technical Journal, 1949, 28(4): 656–715

[29] Sukumar N. Construction of polygonal interpolants: a maximum entropy approach. International Journal for Numerical Methods in Engineering, 2004, 61(12): 2159–2181

[30] Arroyo M, Ortiz M. Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods. International Journal for Numerical Methods in Engineering, 2006, 65(13): 2167–2202

[31] Millán D, Sukumar N, Arroyo M. Cell-based maximum-entropy approximants. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 712–731

[32] Ortiz A, Puso M, Sukumar N. Maximum-entropy meshfree method for compressible and near-incompressible elasticity. Computer Methods in Applied Mechanics and Engineering, 2010, 199(25): 1859–1871

[33] Ortiz A, Puso M, Sukumar N. Maximum-entropy meshfree method for incompressible media problems. Finite Elements in Analysis and Design, 2011, 47(6): 572–585

[34] Quaranta G, Kunnath S K, Sukumar N. Maximum-entropy meshfree method for nonlinear static analysis of planar reinforced concrete structures. Engineering Structures, 2012, 42: 179–189

[35] Ullah Z, Coombs W, Augarde C. An adaptive finite element/meshless coupled method based on local maximum entropy shape functions for linear and nonlinear problems. Computer Methods in Applied Mechanics and Engineering, 2013, 267: 111–132

[36] Amiri F, Anitescu C, Arroyo M, Bordas S P A, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[37] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[38] Amiri F, Millán D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 254–275

[39] Wu C, Young D, Hong H. Adaptive meshless local maximum-entropy finite element method for convection-diffusion problems. Computational Mechanics, 2014, 53(1): 189–200

[40] Kardani O, Nazem M, Kardani M, Sloan S. On the application of the maximum entropy meshfree method for elastoplastic geotechnical analysis. Computers and Geotechnics, 2017, 84: 68–77

[41] Nazem M, Kardani M, Bienen B, Cassidy M. A stable maximum-entropy meshless method for analysis of porous media. Computers and Geotechnics, 2016, 80: 248–260

[42] Navas P, López-Querol S, Yu R C, Li B. Meshfree Methods Applied to Consolidation Problems in Saturated Soils. In: Weinberg K, Pandolfi A, eds. Innovative Numerical Approaches for Multi-Field and Multi-Scale Problems. Springer, 2016, 241–264

[43] Navas P, López-Querol S, Yu R C, Li B. B-bar based algorithm applied to meshfree numerical schemes to solve unconfined seepage problems through porous media. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(6): 962–984

[44] Navas P, Yu R C, López-Querol S, Li B. Dynamic consolidation problems in saturated soils solved through u–w formulation in a LME meshfree framework. Computers and Geotechnics, 2016, 79: 55–72

[45] Zakrzewski N, Nazem M, Sloan S W, Cassidy M On application of the maximum entropy meshless method for large deformation analysis of geotechnical problems. In: Gu Y, Guan H, Sauret E, Saha S, Zhan H, Persky R, eds. Applied Mechanics and Materials. 2016. Trans Tech Publ, 331–335

[46] Lewis R W, Schrefler B A. The finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. John Wiley& Sons, 1998

[47] Jaynes E T. On the rationale of maximum-entropy methods. Proceedings of the IEEE, 1982, 70(9): 939–952

[48] Gull S F, Skilling J. Maximum entropy method in image processing. In: IEE Proceedings F- Communications, Radar and Signal Processing. IET, 1984

[49] Golan A, Judge G G, Miller D. Maximum Entropy Econometrics. John Wiley & Sons, 1996

[50] .Karmeshu J. Entropy Measures, Maximum Entropy Principle and Emerging Applications. Springer Science & Business Media, 2003

[51] Jaynes E T. Information theory and statistical mechanics. Physical Review, 1957, 106(4): 620–630

[52] Gawin D, Baggio P, Schrefler B A. Coupled heat, water and gas flow in deformable porous media. International Journal for Numerical Methods in Fluids, 1995, 20(8–9): 969–987 doi:10.1002/fld.1650200817

[53] Khoei A, Mohammadnejad T. Numerical modeling of multiphase fluid flow in deforming porous media: a comparison between two-and three-phase models for seismic analysis of earth and rockfill dams. Computers and Geotechnics, 2011, 38(2): 142–166

[54] Schrefler B A, Scotta R. A fully coupled dynamic model for two-phase fluid flow in deformable porous media. Computer Methods in Applied Mechanics and Engineering, 2001, 190(24–25): 3223–3246

[55] Brooks R H, Corey A T, 0. Hydraulic properties of porous media and their relation to drainage design. Transactions of the ASAE. American Society of Agricultural Engineers, 1964, 7(1): 26–28

[56] Booker J R, Small J. Finite layer analysis of consolidation. I. International Journal for Numerical and Analytical Methods in Geomechanics, 1982, 6(2): 151–171

[57] Booker J, Small J. A method of computing the consolidation behaviour of layered soils using direct numerical inversion of Laplace transforms. International Journal for Numerical and Analytical Methods in Geomechanics, 1987, 11(4): 363–380

[58] Gibson R, Schiffman R, Pu S. Plane strain and axially symmetric consolidation of a clay layer on a smooth impervious base. Quarterly Journal of Mechanics and Applied Mathematics, 1970, 23(4): 505–520

[59] Aboustit B, Advani S, Lee J. Variational principles and finite element simulations for thermo-elastic consolidation. International Journal for Numerical and Analytical Methods in Geomechanics, 1985, 9(1): 49–69

[60] Liakopoulos A C. Transient Flow Through Unsaturated Porous Media. Dissertation for PhD degree. University of California, Berkeley. 1964

[61] Narasimhan T N, Witherspoon P. Numerical model for saturated-unsaturated flow in deformable porous media: 3. Applications. Water Resources Research, 1978, 14(6): 1017–1034

[62] Schrefler B, Simoni L. A unified approach to the analysis of saturated-unsaturated elastoplastic porous media. Numerical Methods in Geomechanics, 1988, 1: 205–212

[63] Zienkiewicz O, Xie Y M, Schrefler B A, Ledesma A, BicanicN. Static and dynamic behaviour of soils: a rational approach to quantitative solutions. II. Semi-saturated problems. In: Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1990

[64] Schrefler B, Zhan X. A fully coupled model for water flow and airflow in deformable porous media. Water Resources Research, 1993, 29(1): 155–167

[65] Gawin D, Schrefler B A, Galindo M. Thermo-hydro-mechanical analysis of partially saturated porous materials. Engineering Computations, 1996, 13(7): 113–143

[66] .Wang X W, Schrefler B. Fully coupled thermo-hydro-mechanical analysis by an algebraic multigrid method. Engineering Computations, 2003, 20(2): 211–229

[67] Ehlers W, Graf T, Ammann M. Deformation and localization analysis of partially saturated soil. Computer Methods in Applied Mechanics and Engineering, 2004, 193(27): 2885–2910

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature