1 Introduction

This paper introduces a unified formalism for coupled anomalous diffusion processes in porous media as well as a simple open-source implementation in FE commercial code. It was noticed that in soft tissue hydrated with water such as articular cartilage [1] and meniscus [2,3,4,5,6] the interstitial water pressurizes when the tissue is loaded. The fluid pressurization has been hypothesized to be a major factor in the load-support mechanism and in the response to friction of these types of tissues [1]. When the physiological load is applied, this load is shared between the solid skeleton and the interstitial fluid. Theoretical studies have demonstrated that the interstitial fluid (i.e. pore pressure) may sustain up to 90% or more of the total applied load particularly in the initial phases of loading [1, 7, 8].

Experimental studies have revealed that the time-varying response of pore pressure diffusion is anomalous and it is mainly due to variable fluid flow rates which are, in turn, caused by spatially and temporally varying permeability [2]. This is mainly due to the fact that the permeability, hence the rate of fluid flow, often depends on the porosity [9] which varies throughout the tissue [2]. Variations in permeability occur, for example, when the fluid flow impacts the geometry or the micro-structural features such as the configuration of the pores. Experiments on water flow in building materials, sand and zeolite among other minerals, highlighted that the permeability changes during the water flow process as a result of the microstructural rearrangement of grains/pores. Iaffaldano at al. [10] hinted that during compaction of sand, permeability might decrease due to the fact that the fluid carries solid particles which then close some of the pores. Essentially the configuration of the medium, in particular the ratio between closed/open pores, changes during the process. Fluid might be trapped in the medium leading to a slower fluid flow rate. On the contrary, if during the water diffusion process some of the pores open creating conductive microchannels, permeability might increase. Therefore, water can be transported for large distances in a reduced time determining a faster diffusion process. These types of phenomena are well captured by a form of Darcy’s law in which fractional operators are involved in the rate of fluid flow [11, 12].

It is then fundamental to be able to understand and model the time-varying interstitial fluid pressurization, the order of the fractional operators and its role in load bearing of soft hydrated tissues. It is shown that a common mathematical formulation is able to describe both stress-driven anomalous diffusion, biomechanical processes in soft tissues and a number of multiphysics problems. Coupled models of transient flux (temperature, mass, pore pressure) with deformable solids rely on appropriately considering the interaction among separate physical fields. The physics refers to common types of physical processes, e.g., heat transfer (thermo-), pore water or pore pressure movement (hydro/poro-), concentration field (concentro or diffuso/convecto/advecto-), stress and strain (mechano-), dynamics (dyno-), chemical reactions (chemo- or chemico-), electrostatics (electro-), and magnetostatics (magneto-). The mathematics of these phenomena is similar. The differences lay mainly on the coupling terms, and how they are modelled, as they are dictated by the physics of the problem. It is clear that in the case of thermoelasticity and poroelasticity the coupling stress-diffusion in the transport equation is confined to a source/sink term. Different is the case of hydrogen diffusion-mechanics in which the coupling is both in the flux (as there is a term which is related to the gradient of hydrostatic stress) and in a sink term (related to the role of plastic strain).

Here the focus of the paper is to present an unified mathematical and computational framework for mechanics-coupled “anomalous” transport phenomena in porous media and provide the link how the procedure can be applied in a number of multiphysics problems. The anomalous pore pressure diffusion which contain a fractional Darcy’s law is derived and implemented in the finite element framework using the discretised formulation given by Grunwald–Letnikov (GL). The coupled poro-mechanics problem is then, for the first time, implemented in Abaqus, through a UMATHT routine, in order to adopt the coupled temperature-displacement procedure available in Abaqus to effectively solve diffusion-mechanics problems. The same approach has been used when coupling hydrogen diffusion mechanics problems to highlight some material deterioration phenomena. For example the coupled elasto-plastic response with models of diffusional hydrogen transport provided new insights into a number of hydrogen embrittlement mechanisms [13,14,15].

It is also important to know when a given coupled problem can be solved in an uncoupled manner. The strength of coupling is discussed in terms of deriving dimensionless parameters for thermoelasticity/poroelasticity in the same fashion as [16]. Dimensionless parameters are also here derived for the first time for hydrogen diffusion-elasto-plasticity. Summarizing the main novelties of the paper are:

  • Numerical implementation of the anomalous pore pressure diffusion equation in a finite element framework.

  • Details of the Abaqus implementation in UMATHT subroutine and comparison with analytical solutions of a fractional boundary value problems.

  • Derivation of dimensioneless parameters to evaluate the strength of coupling for diffusion mechanics problems

The paper is structured as follows. First the coupled transport-mechanics problem is porous solid is introduced within the context of poroelasticity. The anomalous pore pressure diffusion equation, needed to solve the poroelastic set of equations, is subsequently derived. The strength of coupling in diffusion-mechanics problem is then discussed in Sect. 3. Section 4 focuses on the numerical implementation of the pore pressure diffusion equation, the similarities with heat and hydrogen transport are clearly identified. Numerical results including uncoupled and fully coupled poromechanics such as consolidation problems are presented in Sect. 5.

2 Fully coupled transport and mechanical behaviour in deformable porous solids

In the following we recap the main equations needed to solve coupled pore pressure diffusion with mechanics in the context of the theory of poroelasticity. For the sake of simplicity we limit the discussion to linear elastic solid. In appendix (“Appendix A”) parallelisms with thermo-elasticity and hydrogen diffusion mechanics are reported.

The classical linear model of transient flow and deformation of a homogeneous fully saturated elastic porous medium depends on an appropriate coupling of the fluid pressure and solid stress [17, 18]. A change in applied stress produces a change in fluid pressure or fluid mass and a change in fluid pressure or fluid mass is responsible for a change in the volume of the porous material. The coupling term affects only the hydrostatic part of the stress tensor. Considering that stresses are positive when they are tensile and pressure is positive when it is compressive (i.e. \(tr(\sigma )= -3p\), where p is the pore pressure), the stress tensor is written as follows:

$$\begin{aligned} \varvec{\sigma }= 2G\varvec{\epsilon } + \lambda \text {tr}(\varvec{\epsilon })\mathbf{I} - \alpha p \mathbf{I} \end{aligned}$$
(1)

where \(\lambda = K-\frac{2}{3G}\) is the Lamé constant, GK are the shear and bulk modulus, respectively, p is the pore pressure and \(\alpha \) is the Biot coefficient. It is useful to express the linear poroelastic constitutive equations in terms of strain:

$$\begin{aligned} \varvec{\epsilon }= \frac{\varvec{\sigma }}{2G} - \frac{\nu }{2G(1+\nu )} \text {tr}(\varvec{\sigma })\mathbf{I} - \frac{\alpha }{3K} p\mathbf{I} \end{aligned}$$
(2)

With \(\nu \) being the Poisson’s ratio. It is also important for the following to derive an expression for the volumetric strain or dilatation \(\epsilon _{d}\) by taking the trace of both sides of Eq. 2:

$$\begin{aligned} \epsilon _{d}= \text {tr} (\varvec{\epsilon })= \frac{\text {tr}(\varvec{\sigma })}{2G} - \frac{3\nu }{2G(1+\nu )} \text {tr}(\varvec{\sigma }) - \frac{\alpha }{K} p \end{aligned}$$
(3)

Equation 3 can also be written as follows:

$$\begin{aligned} \epsilon _{d}= \frac{\sigma _{H}}{K} -\frac{\alpha }{K} p \end{aligned}$$
(4)

\(\epsilon _{d}\) being the volumetric/dilatation strain and \(\sigma _{H}=\sigma _{ii}/3\) the hydrostatic component of stress.

Equations 1, A.1 and B.1 in “Appendices A and B” respectively have the same structure. These equations govern coupled stress-diffusion phenomena. The diffusion equations of the pore pressure p, temperature \(\theta \) and hydrogen concentration \(c_{L}\) are derived from the same set of equations (i.e. conservation of pore fluid, energy, mass along with the constitutive equations relating fluxes with gradient of \(p,\theta , c_{L}\) (i.e. Darcy’s, Fourier’s or modified Fick’s laws).

2.1 Anomalous pore pressure diffusion equation

The anomalous pore pressure is controlled by a governing diffusion equation. This equation is found by considering conservation of mass for the pore fluid \(\zeta \), along with a constitutive equation that relates the fluid flux \(\mathbf {j_{p}}\) to the pore pressure gradient through non integer order operators. Mass conservation states that the rate of change of the total mass in a volume \(\Omega \) is equal to the flux \(\mathbf {j_{p}}\) through the surface \(\partial \Omega \) which, in absence of sources or sinks, reads as follows:

$$\begin{aligned} \frac{\partial }{\partial t}\int \limits _{\Omega }{\zeta d \Omega +\int \limits _{\partial \Omega }{(\mathbf {j_{p}}}\cdot {\mathbf {n}}) d \partial \Omega =0} \end{aligned}$$
(5)

where \({\mathbf {n}}\) is the outward normal to \(\partial \Omega \). In the following the linear kinematic/small deformation theory is taken into account. Nonlinear kinematic terms in small deformation regime is often associated with partially saturated conditions (Skempton coefficient B<1), which the theory presented in the paper is not suited for. Considering that the theory presented here is limited to the linear kinematic case and applying the divergence theorem on the surface integral, the following equation is derived:

$$\begin{aligned} \frac{\partial }{\partial t}\int \limits _{\Omega }{(\zeta + \varvec{\nabla } \cdot \mathbf {j_{p}} ) d \Omega =0} \end{aligned}$$
(6)

For Eq. 5 to be true for any region \(\Omega \), the integrand must vanish, i.e.

$$\begin{aligned} \frac{\partial \zeta }{\partial t}+ \nabla \cdot \mathbf {j_{p}} = 0\,. \end{aligned}$$
(7)

Equation 7 represents conservation of mass for the pore fluid. The relationship between the fluid-flux vector \(\mathbf {j_{p}}\) and the pore pressure p is given by a modified version of Darcy’s law involving fractional operators, which in the isotropic case and without inertia’s term takes the form [11, 12]:

$$\begin{aligned} \mathbf {j_{p}}=-{\frac{\varvec{k}}{\mu }} D_{0}^{\beta } \left( {\varvec{\varvec{\nabla } }}p \right) \end{aligned}$$
(8)

where \(\varvec{k}\) is the permeability and \(\mu \) the dynamic viscosity. In the following \({\frac{\varvec{k}}{\mu }}\) will be indicated as \( \varvec{\lambda _{\beta }}\). \( D_{0}^{\beta }\) indicates the Caputo’s fractional time derivative [19]. The Caputo’s fractional derivative of order \(\beta \) of the function \(\varvec{\nabla }p\) is defined as:

$$\begin{aligned} (_{a}D_{t}^{\beta } \varvec{\nabla }p)(t) = \frac{1}{\Gamma (n-\beta )} \int _{a}^{t} \frac{ \varvec{\nabla }p^{n} (\tau )}{(t-\tau )^{\beta +1-n}} d\tau \end{aligned}$$
(9)

The expression is valid for \(n-1<\beta <n\), and \(\Gamma \) is the Euler’s Gamma function. One of the advantage of using the definition of Caputo fractional derivative is that initial conditions in terms of integer order derivatives are obtained. This is important for applications in real physical problems. If \(\alpha \rightarrow n \) the Caputo’s derivative is reduced to classical integer order derivatives; moreover the Caputo fractional derivative of a constant is zero as integer order derivative. Given that in the cases considered in this paper the lower bond of integration \(a=0\), Caputo’s fractional derivative of order \(\beta \) \(_{a}D_{t}^{\beta }\) is indicated as \(D_{0}^{\beta }\).

The fractional derivative method offers the possibility to model with reduced number of parameters all of the anomalous diffusion behaviours by changing the order of the derivative. The drawback is that it is difficult to link the order of the derivatives with microstructural features. In order to understand if the choice of a macroscopic fractional model is sensible,recently approaches based on Bayesian model selection have been developed [20,21,22] and the error in the parameters estimation quantified [20,21,22,23,24].

It is generally assumed that the variation of the pore fluid content \(\zeta \) depends linearly on the hydrostatic component of stress \(\sigma _{H}\) (shear stresses have no influence) and on the pore pressure. Hence \(\zeta \) can be expressed as follows [16]:

$$\begin{aligned} \zeta =\frac{\alpha }{k} \sigma _{H}+ \frac{\alpha }{KB} p \end{aligned}$$
(10)

where B is the Skempton coefficient which essentially gives indication of the increase in fluid pressure of an elastic isotropic porous material subjected to undrained loading. B is related to the drain/undrained bulk moduli: K and \(K_{u}\) and the Biot coefficient \(\alpha \) as follows:

$$\begin{aligned} B=\frac{K_{u}-K}{\alpha K_{u}} \end{aligned}$$
(11)

The value of B is always between 0 and 1. When B is 1, the hydrostatic stress or the dilatation strain is completely transferred into changing pore pressure. When B equals to 0 indicates no change in pore pressure after applying the volumetric stress or strain field. Laboratory studies indicate the value of B depends upon the fluid saturated pore volume of the sample [25]. A value of \(B<1\) indicates a partially saturated media. The theory here considers saturated solid for which B is close to 1. The variation of fluid content with time takes the following form:

$$\begin{aligned} \frac{\partial \zeta }{\partial t}=\frac{\alpha }{k} \frac{\partial \sigma _{H}}{\partial t} + \frac{\alpha }{KB} \frac{\partial p}{\partial t} \end{aligned}$$
(12)

Substituting Eqs. 12 and 8 into Eq. 7 and we obtain the following fractional partial differential equation:

$$\begin{aligned} \frac{\partial p}{\partial t}= \frac{KB}{\alpha } \lambda _{\beta } D_{0}^{\beta } \varvec{\nabla }^{2}p- B\frac{\partial \sigma _{H}}{\partial t} \end{aligned}$$
(13)

In here the following rule of fractional derivatives is implied [19]:

$$\begin{aligned} \nabla \cdot D_{0}^{\beta }\nabla p =D_{0}^{\beta } \varvec{\nabla }^{2}p \end{aligned}$$
(14)

It is important to note that, in case of \(\beta =0\), \(\lambda _{\beta }=\frac{k}{\mu }\) with dimension of \(\frac{[L]^{4}}{[F][T]}\). In the case of \(\beta \ne 0\), \(\lambda _{\beta }\) has dimension of \(\frac{[L]^{4}}{[F][T]^{1-\beta }}\). The pore pressure diffusion equation can also be written in terms of strains. In order to do so we start by combining Eqs. 10, 11 and 4 in order to have an expression relating the variation of fluid content depending on dilatation and pore pressure:

$$\begin{aligned} \zeta =\frac{\alpha ^{2}}{K_{u}-K} p+ \alpha \epsilon _{d} \end{aligned}$$
(15)

Differentiating Eq. 15 with respect to time we obtain:

$$\begin{aligned} \frac{\partial \zeta }{\partial t}= \frac{\alpha ^{2}}{K_{u}-K} \frac{\partial p}{\partial t} + \alpha \frac{\partial \epsilon _{d}}{\partial t} \end{aligned}$$
(16)

Substituting Eqs. 16 and 8 into Eq. 7 we obtain the following fractional partial differential equation in terms of strain:

$$\begin{aligned} \frac{\partial p}{\partial t}= \frac{K_{u}-K}{\alpha ^{2}} \lambda _{\beta } D_{0}^{\beta } \varvec{\nabla }^{2}p- \frac{K_{u}-K}{\alpha }\frac{\partial \epsilon _{d}}{\partial t} \end{aligned}$$
(17)

3 The strength of coupling in diffusion mechanics problems

When dealing with coupled models of transient flux (temperature, mass, pore pressure) with deformable solids it is important to highlight when these problems can be treated in an uncoupled manner. In other words, it is important to quantify when is appropriate to compute stress/strain fields ignoring the transient flux and viceversa. For this purpose in the following dimensionless parameters are identified for the case of poroleasticity, thermoelasticity and hydrogen diffusion-elasto plasticity.

3.1 Poroelastic coupling

Equations 1 and 2 show that the influence of the pore pressure on the volumetric part of the stress or strain tensor is governed by the Biot coefficient \(\alpha \). If \(\alpha =0\) pore pressure would have no influence on the stress/strain tensor. The Biot coefficient is the fluid volume change induced by bulk volume changes in the drained condition. It is shown that, for a large class of porous media, the value of the Biot coefficient \(\alpha \) is significantly greater than zero therefore stresses/strain cannot be computed without computing the pore pressure [16]. Equations 13 and 17 show that the pore pressure is governed by a diffusion-type equation containing an additional coupling term that serves as a source-sink:

$$\begin{aligned} \frac{\partial p}{\partial t}+ \frac{KB}{\alpha }\nabla \cdot \mathbf {j_{p}} +{r}_{p}=0 \end{aligned}$$
(18)

with \({r}_{p}\) being a source term. From Eq. 17:

$$\begin{aligned} Source: {r}_{p}=\frac{K_{u}-K}{\alpha }\frac{\partial \epsilon _{d}}{\partial t}= B K \frac{\partial \epsilon _{d}}{\partial t} \end{aligned}$$
(19)

or from Eq. 13

$$\begin{aligned} Source: {r}_{p}=B\frac{\partial \sigma _{H}}{\partial t} \end{aligned}$$
(20)

These source terms in Eqs. 19 and 20 are the coupling terms relating the pore pressure increment to the increment of the volumetric part of either the stress or strain tensor through the Skempton coefficient B. This effect can be viewed as a transient Skempton-type effect, in the sense that either the change in time of hydrostatic stress or of dilatation strain (i.e. changing in volume) gives rise to a local increase in the pore pressure. Equation 4 relates the dilation strain with hydrostatic stress and pore pressure:

$$\begin{aligned} \frac{\partial \epsilon _{d}}{\partial t}= \frac{1}{K} \frac{\partial \sigma _{H}}{\partial t} -\frac{\alpha }{K} \frac{\partial p}{\partial t} \end{aligned}$$
(21)

Now let us consider two cases:

  1. (a)

    Constant hydrostatic stress: \(\frac{\partial \sigma _{H}}{\partial t}=0\). This reduces Eq. 21 to:

    $$\begin{aligned} \frac{\partial \epsilon _{d}}{\partial t}= -\frac{\alpha }{K} \frac{\partial p}{\partial t} \end{aligned}$$
    (22)

    In the case of constant hydrostatic state of stress, a positive change in pore pressure contributes to lowering the dilation strain rate. Substituting this into Eq. 19 we obtain that the magnitude of source terms becomes:

    $$\begin{aligned} Magnitude Source: \frac{K_{u}-K}{\alpha } \frac{\partial \epsilon _{d}}{\partial t}= \alpha B \frac{\partial p}{\partial t} \end{aligned}$$
    (23)
  2. (b)

    Constant bulk strain:\(\frac{\partial \epsilon _{d}}{\partial t}=0\). This reduces Eq. 21 to:

    $$\begin{aligned} \frac{\partial \sigma _{H}}{\partial t} =\alpha \frac{\partial p}{\partial t} \end{aligned}$$
    (24)

Equation 24 shows that a positive increment of pore pressure causes an increase of the rate of hydrostatic stress. Substituting this into Eq. 20 we obtain that the magnitude of source terms becomes:

$$\begin{aligned} Magnitude Source: B \alpha \frac{\partial p}{\partial t} \end{aligned}$$
(25)

From Eqs. 23 and 25 the strength of coupling is given by \(B \alpha \). If \(B \alpha<<1\) then it is possible to ignore the source terms in the pore pressure diffusion equation written in term of stress Eq. 13 or strain Eq. 17.

It is important to note that substituting Eq. 25 into Eq. 13 we obtain

$$\begin{aligned} \frac{\partial p}{\partial t}= \frac{KB}{\alpha } \lambda _{\beta } D_{0}^{\beta } \varvec{\nabla }^{2}p- \alpha B\frac{\partial p}{\partial t} \end{aligned}$$
(26)

Therefore if \(\alpha B<<1\) then is it possible to ignore mechanical effects. The pore pressure field can then be calculated by solving the standard uncoupled diffusion equation:

$$\begin{aligned} \frac{\partial p}{\partial t}= \frac{KB}{\alpha } \lambda _{\beta } D_{0}^{\beta } \varvec{\nabla }^{2}p \end{aligned}$$
(27)

3.2 Thermoelastic coupling

The temperature diffusion is governed by a diffusion-type equation given below which contains an additional coupling term that serves as a sink (more details in “Appendix C”, i.e. Eqs. C.5 and C.6) :

$$\begin{aligned} \frac{\partial \theta }{\partial t}+ \frac{1}{\rho c_{v}+9K \gamma ^{2}}\nabla \cdot {{\mathbf {j}}}_{\theta }-{r}_{\theta }=0 \end{aligned}$$
(28)

with \({r}_{\theta }\) being a sink term

$$\begin{aligned} Sink: {r}_{\theta }=\frac{3\gamma K \theta }{\rho c_{v}}\frac{\partial \epsilon _{d}}{\partial t} \end{aligned}$$
(29)

and

$$\begin{aligned} Sink: {r}_{\theta }=\frac{3\gamma \theta }{\rho c_{v}+9K \gamma ^{2}\theta }\frac{\partial \sigma _{H}}{\partial t} \end{aligned}$$
(30)

It is clear that the influence of the temperature on volumetric part of the stress or strain tensor is governed by the thermal expansion coefficient \(\gamma \) as stated by Eqs. A.1) and A.2. If \(\gamma =0\) temperature would have no influence on the stress/strain tensor. The coefficient of thermal expansion \(\gamma \) describes how the size of an object changes with a change in temperature. Specifically, it measures the fractional change in size per degree change in temperature at a constant pressure. It is shown that for a wide range of material \(\gamma \) is significantly greater than zero, therefore stresses/strain cannot be computed without computing the temperature [16].

It is possible to relate the dilation strain with hydrostatic stress and temperature (see Eq. A.4 in “Appendix A”):

$$\begin{aligned} \frac{\partial \epsilon _{d}}{\partial t}= \frac{1}{K} \frac{\partial \sigma _{H}}{\partial t} +3 \gamma \frac{\partial \theta }{\partial t} \end{aligned}$$
(31)

As done in the poroelastic coupling section we can consider two cases:

  1. (a)

    Constant hydrostatic stress: \(\frac{\partial \sigma _{H}}{\partial t}=0\). This would mean from Eq. 31 that:

    $$\begin{aligned} \frac{\partial \epsilon _{d}}{\partial t}= 3 \gamma \frac{\partial \theta }{\partial t} \end{aligned}$$
    (32)

    Equation 32 indicates that a positive increment in temperature generates a positive increment in dilation strain, therefore a change in volume of the solid, as it is expected. Substituting this into Eq. 29 we obtain that the magnitude of source terms becomes:

    $$\begin{aligned} \text {Magnitude Sink}: \frac{9K \gamma ^{2} \theta }{\rho C_{v}}\frac{\partial \theta }{\partial t} \end{aligned}$$
    (33)
  2. (b)

    Constant bulk strain:\(\frac{\partial \epsilon _{d}}{\partial t}=0\). This would mean from Eq. 31 that:

    $$\begin{aligned} \frac{\partial \sigma _{H}}{\partial t} = -3K \gamma \frac{\partial \theta }{\partial t} \end{aligned}$$
    (34)

Equation 34 shows that, when the volumetric strain is constant, an increase in temperature is compensated by a decrease in hydrostatic stress. Substituting this into Eq. 30 we obtain that the magnitude of source terms becomes:

$$\begin{aligned} \text {Magnitude Sink}: \frac{9K \gamma ^{2} \theta }{\rho C_{v}+9K \gamma ^{2} \theta } \frac{\partial \theta }{\partial t}=\frac{9K \gamma ^{2} \theta /\rho C_{v}}{1+9K \gamma ^{2} \theta /\rho C_{v}} \frac{\partial \theta }{\partial t} \nonumber \\ \end{aligned}$$
(35)

From Eqs. 33 and 35 the strength of coupling is given by \(\frac{9K \gamma ^{2} \theta }{\rho C_{v}} \). If \(\frac{9K \gamma ^{2} \theta }{\rho C_{v}}<<1\) then it is possible to ignore the sink terms in the temperature diffusion equation.

It is important to note that substituting Eq. 33 into Eq. C.5 we obtain:

$$\begin{aligned} \frac{\partial \theta }{\partial t}= & {} \frac{K_{t}}{\rho c_{v}+9K \gamma ^{2}\theta } \varvec{\nabla }^{2}\theta \nonumber \\&+ \frac{3\gamma \theta }{\rho c_{v}}+\frac{9K \gamma ^{2} \theta /\rho C_{v}}{1+9K \gamma ^{2} \theta /\rho C_{v}} \frac{\partial \theta }{\partial t} \end{aligned}$$
(36)

Therefore, if \(\frac{9K \beta ^{2} \theta }{\rho C_{v}}<<1\) then the temperature field can be calculated by solving the uncoupled diffusion equation for temperature:

$$\begin{aligned} \frac{\partial \theta }{\partial t}= \frac{K_t}{\rho c_{v}} \varvec{\nabla }^{2}\theta \end{aligned}$$
(37)

3.3 Hydrogen diffusion-mechanics coupling

Hydrogen atoms in metals reside either at normal interstitial sites (NILs) or at trapping sites such as dislocations, grain boundaries, carbide/matrix interfaces, microvoids and other defects. We denote the hydrogen concentration in the lattice as \(c_{L}\) (number of H atoms per unit volume) and use \(c_{X}\) to denote the concentration associated with the hydrogen in the traps. The total concentration of hydrogen is given by: \(c_{T}=c_{L}+c_{X}\). The hydrogen diffusion equation reads as follows (“Appendix D”, Eq. D.8):

$$\begin{aligned} \frac{\partial {c}_{L}}{\partial t} + \nabla \cdot \frac{D_{eff}}{D_L} {{\mathbf {j}}}_{c}+{r}_{c}=0 \end{aligned}$$
(38)

with \({r}_{c}\) being a source term, \(c_{L}\) hydrogen concentration at the lattice sites, \({\mathbf {j}}_{c}=-\frac{{{D}_{L}}{{c}_{L}}}{RT}\nabla \mu \) hydrogen flux. \(\mu =\mu _{0}+ RT\ln {c}_{L} -\sigma _{H} V_H\) is the chemical potential. The source terms is expressed as below:

$$\begin{aligned} Source: {r}_{c}=\frac{D_{eff}}{D_L} \frac{\partial {c}_{X}}{\partial {N}_X} \frac{d{N}_X}{d\varepsilon ^p} \frac{d \varepsilon ^{p}}{d t} \end{aligned}$$
(39)

where \(N_{X}\) is the density of traps which increases with increasing plastic strain.

Hydrogen concentration influences both (a) the volumetric part of the stress/strain tensor and (b) plastic deformation.

  1. (a)

    The influence of hydrogen concentration on the volumetric part of the stress or strain tensor is ruled by Eqs. B.1 and B.2 reported in “Appendix B”. It can be noted that if \(V_{H} c_{L}=0\) hydrogen concentration would have no influence on the stress/strain tensor. The partial molar volume of hydrogen \(V_{H}\) is related to the occupancy of available hydrogen sites. If there are no sites to be occupied by hydrogen atoms then hydrogen concentration has no influence on the stress tensor. Summarizing, the dependence on hydrostatic stress can be neglected if \(V_{H} c_{L}<<1\). \(V_{H} c_{L}\) is an adimensional parameter that can be used to quantify the coupling with the elastic stress.

  2. (b)

    The influence of plastic strain in H-diffusion process can be neglected if the source term in Eq. 39 is negligible. This happens if \(\frac{\partial {c}_{X}}{\partial {N}_X}<<1\) when traps binding energy is high as \(\frac{\partial {c}_{X}}{\partial {N}_X}= \frac{ K_c c_{L}}{\beta N_{L}+K_c c_{L}}\) and \(K_c=\exp (-\frac{W_{\beta }}{RT})\) with \(W_{\beta }\) being the binding energy.

So if \(-\frac{RT}{W_{\beta }}<<1\) then dependence of plastic strain on H-transport can be neglected. The adimensional parameter that can be used to quantify the coupling with the plastic strain is: \(-\frac{RT}{W_{\beta }}\).

If both terms: \(V_{H} c_{L}\) and \(-\frac{RT}{W_{\beta }}\) are close to zero then the hydrogen diffusion equation to be solved is as follow:

$$\begin{aligned} \frac{\partial {c}_{L}}{\partial t} (1+ \frac{\partial {c}_{X}}{\partial {c}_{L}}) - \nabla \cdot (D_{L} \nabla {c}_{L}) =0 \end{aligned}$$
(40)

4 Analogy between coupled transport equations and numerical implementation

The form of the transport equations in Eqs. 18, 28 and 38 is similar and they can be summarized here below in a generic form:

$$\begin{aligned} \frac{\partial X}{\partial t} + Y \nabla \cdot J + R=0 \end{aligned}$$
(41)

where X is the degree of freedom which is \(\theta \) for the the heat diffusion, \({c_{L}}\) for hydrogen transport and p for pore pressure diffusion equations. The similarities of these equations are summarized in Table 1.

Table 1 Analogy of variables between heat transfer, mass diffusion and pore pressure diffusion

The transport equation in Eq. 41 in the frame of finite element method (FEM) is solved step by step with numerical integration following an implicit scheme (Newton-Raphson algorithm). Considering that:

$$\begin{aligned} \frac{\partial X}{\partial t}= {\dot{X}}= \frac{X_{t+\Delta _{t}} - X_{t}}{\Delta _{t}} \end{aligned}$$
(42)

Equation 41 reads as follow:

$$\begin{aligned} \int \limits _{\Omega }{ {\dot{X}} }d\Omega =-\int \limits _{\Gamma }{Y({{\mathbf {J}}\cdot {\mathbf {n}}})} d\Gamma + \int \limits _{\Omega }{R}d\Omega =0 \end{aligned}$$
(43)

Applying the divergence theorem and writing the weak form Eq. 43 becomes:

$$\begin{aligned} \int \limits _{\Omega }{\delta X{\dot{X}}} d\Omega = - \int \limits _{\Omega } Y\delta X (\nabla \cdot {{\mathbf {J}}}) d\Omega + \int \limits _{\Omega } {\delta {X} R} d\Omega =0 \end{aligned}$$
(44)

\(\delta X\) is an arbitrary variational field. Considering that:

$$\begin{aligned} \int \limits _{\Omega } Y \delta (\nabla \cdot X {{\mathbf {J}}}) d\Omega = \int \limits _{\Gamma }{ Y \delta (X{{\mathbf {J}}\cdot {\mathbf {n}}})} d\Gamma \end{aligned}$$
(45)

and

$$\begin{aligned} q=- {\mathbf {J}}\cdot {\mathbf {n}} \end{aligned}$$
(46)

Equation 44 becomes:

$$\begin{aligned} \int \limits _{\Omega }{\delta X{\dot{X}}} d\Omega= & {} \int \limits _{\Omega } Y \delta (\nabla X \cdot {{\mathbf {J}}}) d\Omega + \int \limits _{\Gamma }{Y \delta X (q)} d\Gamma \nonumber \\&\quad + \int \limits _{\Omega } {Y \delta X R} d\Omega =0 \end{aligned}$$
(47)

Substituting Eq. 42 into Eq. 47 the following relation is derived:

$$\begin{aligned}&\int \limits _{\Omega }{\delta X \frac{X_{t+\Delta _{t}} - X_{t}}{\Delta _{t}}} d\Omega = \int \limits _{\Omega } Y \delta (\nabla X \cdot {{\mathbf {J}}}) d\Omega \nonumber \\&\quad + \int \limits _{\Gamma }{Y \delta X (q)} d\Gamma + \int \limits _{\Omega } {\delta X R} d\Omega =0 \end{aligned}$$
(48)

Now we indicate with U: :

$$\begin{aligned} U=\frac{\partial X}{\partial \Delta _{t}}= U (X, t, \nabla X, S^{i}) \end{aligned}$$
(49)

\(S^{i}\) state variables.

$$\begin{aligned} J= J (X, t, \nabla X, S^{i}) \end{aligned}$$
(50)

and source/sink

$$\begin{aligned} R=R (X, t, \nabla X, S^{i}) \end{aligned}$$
(51)

4.1 Numerical Implementation of the fractional pore pressure diffusion equation in the UMATHT subroutine

The pore pressure diffusion equation (Eq. 18) can be solved by substituting \(X=p, J=j_{p}, R=r_{p}, U=U_{p}, Y=\frac{KB}{\alpha }\) in Eq. 48 to obtain:

$$\begin{aligned}&\int \limits _{\Omega }{U_{p} (p_{t+ \Delta _{t}}- p_{t}) } d\Omega = \int \limits _{\Omega } {\delta \nabla p \cdot {{\mathbf {J}}}_{p}} d\Omega + \int \limits _{\Gamma }{\delta p (q)} d\Gamma \nonumber \\&\quad + \int \limits _{\Omega } {\delta _{p} r_{p}} d\Omega =0 \end{aligned}$$
(52)

A list of the terms to be coded in UMATHT routine are given below:

\(U_{p}\):

\(U_{p}\):

$$\begin{aligned} U_{p}(t+\Delta _{t})=U_{p}(t)+\frac{\partial U_{p}}{\partial p}dp+\frac{\partial r_{p}}{\partial \sigma _{H}} d\sigma _{H} \end{aligned}$$
(53)

where:

$$\begin{aligned}&\frac{\partial U_{p} }{\partial p} =1 \end{aligned}$$
(54)
$$\begin{aligned}&\frac{\partial U_{p} }{\partial \nabla p} =0 \end{aligned}$$
(55)
  • Source:

    $$\begin{aligned} \frac{\partial r_{p}}{\partial \sigma _{H}} d\sigma _{H} \end{aligned}$$
    (56)

    where:

    $$\begin{aligned} \frac{\partial r_{p}}{\partial \sigma _{H}} =-B \end{aligned}$$
    (57)
  • Flux: The flux term in Eq. 8 contains a fractional derivative. The fractional derivatives must be discretized. To this purpose the Grunwald–Letnikov (GL) fractional derivative is used [19]:

    $$\begin{aligned}&\left( _{0}^{GL}D_{t}^{\beta } \nabla p \right) (t)=\left( _{0}^{GL}D_{t}^{\beta } \nabla p \right) (k\Delta _{t})\nonumber \\&\quad =\lim _{\Delta _{t}\rightarrow 0}\,\Delta t^{-\beta }\sum _{j=1}^{k+1}\lambda _{j}^{(\rho )} \nabla p^{(k-j+1)} \end{aligned}$$
    (58a)
    $$\begin{aligned}&\lambda _{j+1}^{(\beta )}=\frac{j-1-\beta }{j}\lambda _{j};\;\;\;\;\;\;\;\lambda _1=1 \end{aligned}$$
    (58b)

    where \(j_{p}^{(k-j+1)}=j_{p}[(k-j+1)\Delta _{t}]\). For sufficiently small \(\Delta _{t}\) the GL fractional derivative coalesces with the Caputo’s fractional derivative.

    $$\begin{aligned}&j_{p}^{k} =-\frac{KB}{\alpha }\lambda _{\beta } \Delta _{t}^{-\beta } \sum _{j=1}^{k} c_{j} \nabla p^{k-j+1} \end{aligned}$$
    (59)
    $$\begin{aligned}&\frac{\partial j_{p} }{\partial \nabla p} =-\frac{KB}{\alpha }\lambda _{\beta } \Delta _{t}^{-\beta } \end{aligned}$$
    (60)
  • Jacobian:

    $$\begin{aligned}&\frac{1}{\Delta _{t}} \int \limits _{\Omega }{\delta p \frac{\partial U_{p} }{\partial p} \cdot dp} d\Omega +\frac{1}{\Delta _{t}} \int \limits _{\Omega } {\delta p \frac{\partial U_{p} }{\partial \nabla p} \cdot dp} d\Omega \nonumber \\&\quad - \int \limits _{\Omega } {\delta \nabla p \cdot \frac{\partial j_{p} }{\partial p} dp} d\Omega \nonumber \\&- \int \limits _{\Omega } {\delta \nabla p \cdot \frac{\partial j_{p} }{\partial \nabla p} d \nabla p} d\Omega - \int \limits _{\Omega } {\delta p \frac{\partial r_{p}}{\partial p} \cdot dp} d\Omega \nonumber \\&\quad - \int \limits _{\Gamma }{\delta p \frac{\partial q}{\partial p}\cdot dp } d\Gamma =0 \end{aligned}$$
    (61)

Usually FE codes that use an implicit Newton-Raphson integration scheme allow the time increment to be determined automatically to optimize the run time. The GL formula for evaluating of the fractional derivatives has been derived assuming a constant increment (i.e. the time) and, to the best of our knowledge, a corresponding formulation for a variable increment is not available in the literature; furthermore, the automatic time increment requires the definition of a tolerance criterion, that is difficult to define without knowledge of the elastic and inelastic parts of the strain. For these two reasons we have currently limited ourselves to using this model with a fixed time increment. In order to evaluate the GL derivative the history of strain at each Gauss Point must be stored leading possibly to a considerable amount of memory when analysing large FE models. A number of strategies to overcome this problem have been explored and are discussed in [26].

5 Numerical results

5.1 Uncoupled problem and comparison with analytical solutions

The numerical solution of a simple 1D pore pressure diffusion problem can compare with the analytical solution. A bar made of porous material of length l = 50 mm subjected to an initial pore pressure \(p=0\) Pa throughout the bar. At time t = 0, the pressure at one end is increased to \(p=100\) Pa. The governing fractional pore pressure diffusion equation Eq. 13 (\(\sigma _{H}=0\)) and the boundary conditions are written below:

$$\begin{aligned}&\frac{\partial p}{\partial t}= \frac{KB}{\alpha } \lambda _{\beta } D_{0}^{\beta } \nabla ^{2}p \end{aligned}$$
(62)
$$\begin{aligned}&p(x,0)=0; 0<x<l \end{aligned}$$
(63)
$$\begin{aligned}&p(0,t)=p_{0}=0; p(l,t)=p_{l}=100Pa \end{aligned}$$
(64)

The analytical solution was also found (see “Appendix E”) and it reads as follow:

$$\begin{aligned}&p(x,t)= p_{0} + \frac{x}{l} (p_{l}- p_{0})+ \sum \limits _{n=1}^\infty E_{1-\beta ,1}\nonumber \\&\quad \left[ \left( -{\frac{n^{2}\pi ^{2} {\bar{\lambda }} t^{1-\beta }}{l^{2}}}\right) \right] c_{n}\sin \frac{n \pi x}{l} \end{aligned}$$
(65)

and

$$\begin{aligned} c_{n}= \frac{2}{n \pi } [-1^{n} (p_{l}- p_{0})] \end{aligned}$$
(66)

Note that in Eq. 65 in the case of \(\beta =0\) the Mittag Leffler function \(E_{1-\beta ,1}\) become \(E_{1,1}=\exp \), the exponential function. In this case the solution in Eq. 65 is identical to the classical Fick’s diffusion solution [13]:

$$\begin{aligned}&p(x,t)= p_{0} + \frac{x}{l} (p_{l}- p_{0})+ \sum \limits _{n=1}^\infty \exp \nonumber \\&\quad \left[ \left( -{\frac{n^{2}\pi ^{2} {\bar{\lambda }} t^{1-\beta }}{l^{2}}}\right) \right] c_{n}\sin \frac{n \pi x}{l} \end{aligned}$$
(67)

where \({\bar{\lambda }}=\frac{KB}{\alpha } \lambda _{\beta }\). The material parameters used in the simulations are: \(K=13\times 10^{9}\, \hbox {Pa}, B=0.88, \alpha =0.65, \lambda _{\beta }=8.33\times 10^{-8} \frac{\hbox {m}^{2}}{\hbox {Pa s}}\). The transient solution is affected by the value of \(\beta \). This phenomenon can be appreciated by considering the characteristic time for diffusion i.e. time necessary for the pore pressure to diffuse along the bar. This is calculated from the third term in Eq. 65 being equal to zero:

\(\sum \nolimits _{n=1}^\infty E_{1-\beta ,1} [\left( -{\frac{n^{2}\pi ^{2} {\bar{\lambda }} t^{1-\beta }}{l^{2}}}\right) ] c_{n}\sin \frac{n \pi x}{l}=0\). This equality is reached when the characteristic time \({\bar{t}}\) is equal to \({\bar{t}}=(\frac{l^{2}}{{\bar{\lambda }}})^{\frac{1}{1-\beta }}\).

The steady state solution reads: \(p(x,t)= p_{0} + \frac{x}{l} (p_{l}- p_{0})\).

Figure 1 shows the evolution of the characteristic time \({\bar{t}}=(\frac{l^{2}}{{\bar{\lambda }}})^{\frac{1}{1-\beta }}\) with \(\beta \). In the case of \(\beta =0\) the characteristic time of the standard diffusion case is recovered, i.e. \({\bar{t}}=(\frac{l^{2}}{{\bar{\lambda }}})\). Increasing \(\beta \) (with \(0<\beta <1\)) the characteristic time decreases which means that the transient diffusion phenomenon is faster, therefore is it important to consider the time step to adopt when solving the mathematical problem. In order to capture the transient phenomenon when increasing the value of \(\beta \), the time step will need to be reduced.

Fig. 1
figure 1

Characteristic time for diffusion, i.e. time necessary for pore pressure to diffuse along the bar

Figure 2 shows the evolution of the pore pressure along the bar for \(\beta =0, 0.1, 0.3, 0.5\). The Abaqus simulations were conducted using quadratic continuum coupled temperature-displacements hybrid element CPE8HT. It can be seen that the full transient solutions for the case of \(\beta =0, 0.1, 0.3, 0.5\) at different times are in perfect agreement. The input files and the UMATHT subroutine are attached as electronic supplement material.

Fig. 2
figure 2

Transient evolution of pore pressure: comparison between analytic solutions (continuous line) and simulations (dotted lines). a \(\beta =0\), time step \(dt=1.7\times 10^{-9}\), total time \(t=1.7\times 10^{-7}s\), b \(\beta =0.1\), time step \(dt=1.7\times 10^{-10}\), total time \(t=1.7\times 10^{-8}s\), c \(\beta =0.3\), time step \(dt=1.7\times 10^{-11}\), total time \(t=1.7\times 10^{-9}s\) and d \(\beta =0.5\), time step \(dt=1.7\times 10^{-14}\), total time \(t=1.7\times 10^{-12}s\)

5.2 Fractional consolidation: comparison with analytical solutions

We solve here a fully coupled poromechanics problem. Let us consider a confined consolidation experiment in which we apply a compressive stress on a cylinder made of porous material saturated with water as pictured in Fig. 3. The cylinder is insulated except at the base through which water can flow out. During the deformation process the applied load is sustained by both the pore pressure and the solid skeleton. As for the Terzaghi’s solution, the applied load is firstly sustained by the pore pressure. As fluid flows out of the sample the pore pressure decreases and the solid skeleton starts deforming. We are interested here to model the transient pore pressure diffusion phenomenon and to understand the effect of the order of the fractional derivative \(\beta \).

The test is modelled through the 1D uniaxial strain poroelastic problem. Considering that there is variation of field quantities only in z-direction, the equilibrium equation is given by \(\frac{\partial \sigma _{zz}}{\partial z}=0\), recalling Eq. 1 adapted for the 1D case in which the only non zero component of strain is \(\epsilon _{zz}\), we obtain:

$$\begin{aligned} \left( K+\frac{4G}{3}\right) \frac{\partial \epsilon _{zz}}{\partial z} - \alpha \frac{\partial p}{\partial z}=0 \end{aligned}$$
(68)

Equation 68 is then coupled with the pore pressure diffusion equation in Eq. 17 adapted for the 1D uniaxial strain version below:

$$\begin{aligned} \frac{\partial p}{\partial t}= \frac{K_{u}-K}{\alpha ^{2}} \lambda _{\beta } D_{0}^{\beta } \varvec{\nabla }^{2}p- \frac{K_{u}-K}{\alpha }\frac{\partial \epsilon _{zz}}{\partial t} \end{aligned}$$
(69)

Combining Eqs. 68 and 69 we obtain the following pore pressure diffusion equation:

$$\begin{aligned} \frac{\partial p}{\partial t}= {\bar{\lambda }} D_{0}^{\beta }\frac{\partial ^ {2}p}{\partial z^{2}} \end{aligned}$$
(70)

where \({\bar{\lambda }}= \lambda _{\beta }\frac{\left( 4G+3K\right) \left( K_{u}-K\right) }{\alpha ^{2}\left( 4G+3K_{u}\right) }\). The boundary conditions are those shown in Fig. 3:

$$\begin{aligned}&\frac{\partial p}{\partial z}_{z=0}=0, \end{aligned}$$
(71)
$$\begin{aligned}&p(z=h,t)=0, \end{aligned}$$
(72)
$$\begin{aligned}&u(z=h,t)=0. \end{aligned}$$
(73)
Fig. 3
figure 3

a Schematic representation of the consolidation test in which a compressive force is applied to a saturated sample of height z = 3 mm. All boundaries are impermeable except at the base through which the fluid flows out. b Detail of the discretized domain

A constant compressive stress in the z-direction is applied to the cylinder at \(z=0\):

$$\begin{aligned} \sigma _{zz}(0,t)=- P_{A} \end{aligned}$$
(74)

where \(-P_{A}\) is the applied compressive stress. The initial pore pressure is derived for undrained conditions i.e.:

$$\begin{aligned} p(z,0)= P_{A}\frac{3(K_{u}-K)}{\alpha \left( 4G+3K_{u}\right) } \end{aligned}$$
(75)

The analytical solution in terms of pore pressure reads as follows:

$$\begin{aligned} p(z,t)= P_{A}\gamma \sum \limits _{n=1,3}^\infty E_{1-\beta ,1} \left( -{\frac{n^{2}\pi ^{2} {\bar{\lambda }} t^{1-\beta }}{4h^{2}}}\right) c_{n}\cos \frac{n \pi z}{2h} \nonumber \\ \end{aligned}$$
(76)

where:

$$\begin{aligned} \gamma&= \frac{3\left( K_{u}-K\right) }{\alpha \left( 4G+3K_{u}\right) } \end{aligned}$$
(77)
$$\begin{aligned} c_{n}&= \frac{4}{n \pi } \left( -h\right) ^{\frac{n-1}{2}} \end{aligned}$$
(78)

where \(E_{1-\beta ,1}\) is the Mittag–Leffler function. In the case of \(\beta =0\) the solution identical to the classical Terzaghi’s solution in which \(E_{1-\beta ,1}=\exp \). As done in the previous section, it is interesting to observe that the characteristic time for diffusion for this problem is equal to \({\bar{t}}=(\frac{4h^{2}}{{\bar{\lambda }}})^{\frac{1}{1-\beta }}\). This is calculated as the time needed for the transient diffusion process to finish i.e. the third term in Eq. 76\(\sum \nolimits _{n=1,3}^\infty E_{1-\beta ,1} \left( -{\frac{n^{2}\pi ^{2} {\bar{\lambda }} t^{1-\beta }}{4h^{2}}}\right) =0\).

It is possible to understand the role of the order of the fractional derivative \(\beta \) in the transient solution by observing Fig. 4 in which the variation of the characteristic time \({\bar{t}}\) with \(\beta \) is shown. It can be noted that, as in the pure diffusion case, the transient phenomenon is faster when the order of the fractional derivative \(\beta \) increases, i.e. the characteristic time decreases. However, compared to the pure diffusion case in Fig. 1, the rate of the decrease in the duration of the transient phenomenon with increasing \(\beta \) is slower in the consolidation problem. Figure 5 shows the evolution of the pore pressure along the cylinder for \(\beta =0, 0.1, 0.5\). The Abaqus simulations were conducted using quadratic continuum coupled temperature-displacements hybrid element CPE8HT. The material parameters used in the simulations are summarized in Table 2. In this case we keep the time step and the total time constant (\(dt=4.8\times 10^{-6}\), total time \(t=4.8\times 10^{-4}\)) for the three simulations with the value of \(\beta =0, 0.1, 0.5\). Figure 5a–c show that the rate of pore pressure diffusion increases with increasing order of fractional derivative \(\beta \). This also implies that the rate at which the fluid flows out of the sample during the consolidation problem is higher when the value of \(\beta \) increases.

Fig. 4
figure 4

Characteristic time for the pore pressure diffusion in the case of a consolidation problem

Table 2 Material parameters used to solve the consolidation problem
Fig. 5
figure 5

Transient evolution of pore pressure:comparison between analytic solutions (continuous line) and simulations (dotted lines). a \(\beta =0\), b \(\beta =0.1\), c \(\beta =0.5\), time step \(dt=4.8\times 10^{-6}\), total time \(t=4.8\times 10^{-4}s\)

6 Conclusions

This paper presents a unified formalism for coupled anomalous diffusion processes in porous media as well as a simple open-source implementation within FE software. Here the anomalous pore pressure diffusion, which contains a Darcy’s law involving non-integer order operators, is derived and implemented in the finite element framework using the discretised formulation given by Grunwald–Letnikov (GL). The coupled poro-mechanics problem is then, for the first time, implemented in a UMATHT routine in the commercial software ABAQUS. The similarity of the formalism between coupled temperature-displacement and pore pressure- displacement procedures is used. The unified framework presented here is adaptable for other multiphysics problems ruled by the same set of PDEs such as thermoelasticity and hydrogen diffusion mechanics. The differences among multiphysics problems lay on how the coupling terms are modelled as these are associated with the specific physics of the phenomena. Here we have summarized the fundamental equations for coupling, implementing and solving diffusion equations involved in fractional poroelasticity/thermoelasticty and hydrogen diffusion-mechanics highlighting how the coupling between stress-strain and transport differs in these three cases and the strength of the coupling terms. We have derived dimensionless parameters that allow to understand when is appropriate to simplify and uncouple the problem, i.e. solving the diffusion equation without consider the stress state and viceversa. Numerical simulations of uncoupled and coupled poromechanics problems are successfully compared with analytical solutions of fractional boundary value problems whose details are also given in the appendix.