A finite-strain hyperviscoplastic model and undrained triaxial tests of peat

TL;DR

This paper develops a finite-strain hyperviscoplastic model for peat, incorporating rate-dependent and hysteresis effects, and validates it against undrained triaxial tests at multiple strain rates.

Contribution

It introduces a thermodynamically consistent hyperviscoplastic model for peat and demonstrates its effectiveness through experimental validation under undrained conditions.

Findings

Model accurately simulates undrained triaxial tests at various strain rates.

Experimental data supports the model's ability to capture hysteresis and rate effects.

Model parameters related to solid matrix properties were successfully fitted.

Abstract

This paper presents a finite-strain hyperviscoplastic constitutive model within a thermodynamically consistent framework for peat which was categorised as a material with both rate-dependent and thermodynamic equilibrium hysteresis based on the data reported in the literature. The model was implemented numerically using implicit time integration and verified against analytical solutions under simplified conditions. Experimental studies on the undrained relaxation and loading-unloading-reloading behaviour of an undisturbed fibrous peat were carried out to define the thermodynamic equilibrium state during deviatoric loading as a prerequisite for further modelling, to fit particularly those model parameters related to solid matrix properties, and to validate the proposed model under undrained conditions. This validation performed by comparison to experimental results showed that the hyperviscoplastic model could simulate undrained triaxial compression tests carried out at five different strain rates with loading/unloading relaxation steps.

Tables

Table 1. Table 1: Selective properties of the tested peat material [ 63 ] .

Natural water content per dry mass [%]	Specific gravity of solid constituents	Loss on ignition [%]	Fibre content [%]	pH
590	1.42	98.6	63.5	3.8

Table 2. Table 2: Parameter sensitivity test of the hyperelastoplastic model.

	$C_1$ [kPa]	$D_2$, $D_{2p}$ [kPa]	$\alpha$	$C_{1p}$ [kPa]	${\alpha }_p$	$c_p$ [kPa-1]
Set $C_1$	9.0$\pm$5.0	500.0	0.0	50.0	0.0	0.1
Set $D_2/D_{2p}$	9.0	1000.0$\pm$ 500.0	0.0	50.0	0.0	0.1
Set $\alpha$	9.0	500.0	0.5$\pm$0.5	50.0	0.0	0.1
Set $C_{1p}$	9.0	500.0	0.0	50.0$\pm$20.0	0.0	0.1
Set ${\alpha }_p$	9.0	500.0	0.0	50.0	100.0$\pm$100.0	0.1
Set $c_p$	9.0	500.0	0.0	50.0	0.0	0.1$\pm$0.02

Table 3. Table 3: Material parameters of the hyperelastoplastic model for the tested peat in undrained condition.

$C_1$ [kPa]	$D_2$ [kPa]	$\alpha$	$C_{1p}$ [kPa]	$D_{2p}$ [kPa]	${\alpha }_p$	$c_p$ [kPa-1]
9.0	500.0	0.0	50.0	500.0	0.0	0.1

Table 4. Table 4: Material parameters of the hyperviscoelastic model for the undrained peat.

Parameter	$C_{1v}^i$ [kPa]	$D_{2v}^i$ [kPa]	${\alpha }_v^i$	${\eta }_v^i$ [kPa d]
$i=1$	8.0	500.0	0.0	9.0
$i=2$	40.0	500.0	0.0	0.35

Equations

$$\mathbf{S}={\mathbf{S}}_{\mathrm{eq}}+{\mathbf{S}}_{\mathrm{ov}}$$

$$\sigma =\underset{\text{equilibrium state}}{\underbrace{{\sigma }_{eq0}+{\sigma }_{\mathrm{p}}}}+\sum _{i=1}^2{\sigma }_{\mathrm{ov}}^i$$

$$\epsilon ={\epsilon }_{\mathrm{ep}}+{\epsilon }_{\mathrm{p}}={\epsilon }_{\mathrm{ev}}^i+{\epsilon }_{\mathrm{v}}^i$$

$$\mathbf{S}={\mathbf{S}}_{\mathrm{eq}}+{\mathbf{S}}_{\mathrm{p}}+\sum _{i=1}^2{\mathbf{S}}_{\mathrm{ov}}^i$$

$$\mathbf{F}={\mathbf{F}}_{\mathrm{ep}}{\mathbf{F}}_{\mathrm{p}}={\mathbf{F}}_{\mathrm{ev}}^i{\mathbf{F}}_{\mathrm{v}}^i$$

$$\mathbf{E}=\frac{1}{2}\left({\mathbf{F}}^{\mathrm{T}}\mathbf{F}-\mathbf{I}\right)=\frac{1}{2}\left(\mathbf{C}-\mathbf{I}\right)$$

$$\mathbf{C}={\mathbf{F}}^{\mathrm{T}}\mathbf{F}$$

$$ϵ={\mathbf{F}}_{\mathrm{p}}^{-\mathrm{T}}\mathbf{E}{\mathbf{F}}_{\mathrm{p}}^{-1}=\frac{1}{2}\left({\mathbf{F}}_{\mathrm{ep}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{ep}}-{\mathbf{F}}_{\mathrm{p}}^{-\mathrm{T}}{\mathbf{F}}_{\mathrm{p}}^{-1}\right)={ϵ}_{\mathrm{ep}}+{ϵ}_{\mathrm{p}}$$

$$ϵ={\mathbf{F}}_{\mathrm{v}}^{i-\mathrm{T}}\mathbf{E}{\mathbf{F}}_{\mathrm{v}}^{i-1}=\frac{1}{2}\left({\mathbf{F}}_{\mathrm{ev}}^{i\mathrm{T}}{\mathbf{F}}_{\mathrm{ev}}^i-{\mathbf{F}}_{\mathrm{v}}^{i-\mathrm{T}}{\mathbf{F}}_{\mathrm{v}}^{i-1}\right)={ϵ}_{\mathrm{ev}}^i+{ϵ}_{\mathrm{v}}^i$$

$${ϵ}_{\mathrm{ep}}=\frac{1}{2}\left({\mathbf{F}}_{\mathrm{ep}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{ep}}-\mathbf{I}\right)$$

$${ϵ}_{\mathrm{ev}}^i=\frac{1}{2}\left({\mathbf{F}}_{\mathrm{ev}}^{i\mathrm{T}}{\mathbf{F}}_{\mathrm{ev}}^i-\mathbf{I}\right)$$

$$\stackrel{\Delta }{ϵ}={\mathbf{F}}_{\mathrm{p}}^{-\mathrm{T}}\dot{\mathbf{E}}{\mathbf{F}}_{\mathrm{p}}^{-1}=\dot{ϵ}+ϵ{\mathbf{l}}_{\mathrm{p}}+{\mathbf{l}}_{\mathrm{p}}^{\mathrm{T}}ϵ$$

$$\stackrel{\Delta }{ϵ}={\mathbf{F}}_{\mathrm{v}}^{i-\mathrm{T}}\dot{\mathbf{E}}{\mathbf{F}}_{\mathrm{v}}^{i-1}={\dot{ϵ}}^i+{ϵ}^i{\mathbf{l}}_{\mathrm{v}}^i+{\mathbf{l}}_{\mathrm{v}}^{i\mathrm{T}}{ϵ}^i$$

$${\mathbf{l}}_{\mathrm{p}}$$

$$\stackrel{\Delta }{ϵ}={\stackrel{\Delta }{ϵ}}_{\mathrm{ep}}+{\stackrel{\Delta }{ϵ}}_{\mathrm{p}}$$

$${\stackrel{\Delta }{ϵ}}_{\mathrm{ep}}={\dot{ϵ}}_{\mathrm{ep}}+{ϵ}_{\mathrm{ep}}{\mathbf{l}}_{\mathrm{p}}+{\mathbf{l}}_{\mathrm{p}}^{\mathrm{T}}{ϵ}_{\mathrm{ep}}$$

$${\stackrel{\Delta }{ϵ}}_{\mathrm{p}}={\dot{ϵ}}_{\mathrm{p}}+{ϵ}_{\mathrm{p}}{\mathbf{l}}_{\mathrm{p}}+{\mathbf{l}}_{\mathrm{p}}^{\mathrm{T}}{ϵ}_{\mathrm{p}}=\frac{1}{2}\left({\mathbf{l}}_{\mathrm{p}}+{\mathbf{l}}_{\mathrm{p}}^{\mathrm{T}}\right)$$

$${\mathbf{S}}_{\mathrm{p}}:\dot{\mathbf{E}}={\tau }_{\mathrm{p}}:\stackrel{\Delta }{ϵ}={\tau }_{\mathrm{p}}:\left({\stackrel{\Delta }{ϵ}}_{\mathrm{ep}}+{\stackrel{\Delta }{ϵ}}_{\mathrm{p}}\right)$$

$${\mathbf{S}}_{\mathrm{ov}}^i:{\dot{\mathbf{E}}}^i={\tau }_{\mathrm{ov}}^i:\stackrel{\Delta }{ϵ}={\tau }_{\mathrm{ov}}^i:\left(\stackrel{\Delta }{ϵ}{}_{\mathrm{ev}}^i+\stackrel{\Delta }{ϵ}{}_{\mathrm{v}}^i\right)$$

$${\tau }_{\mathrm{p}}={\mathbf{F}}_{\mathrm{p}}{\mathbf{S}}_{\mathrm{p}}{\mathbf{F}}_{\mathrm{p}}^{\mathrm{T}}$$

$${\tau }_{\mathrm{ov}}={\mathbf{F}}_{\mathrm{v}}^i{\mathbf{S}}_{\mathrm{ov}}^i{\mathbf{F}}_{\mathrm{v}}^{i\mathrm{T}}$$

$$\bar{\psi }\left(\mathbf{E},{ϵ}_{\mathrm{ep}},{ϵ}_{\mathrm{ev}}\right)={\bar{\psi }}_0\left(\mathbf{E}\right)+{\bar{\psi }}_{\mathrm{p}}\left({ϵ}_{\mathrm{ep}}\right)+\sum _{i=1}^2{\bar{\psi }}_{\mathrm{ev}}^i\left({ϵ}_{\mathrm{ev}}\right)$$

$$\mathbf{S}:\dot{\mathbf{E}}-{\rho }_0\dot{\bar{\psi }}\ge 0$$

$$\left({\mathbf{S}}_{\mathrm{eq}}-{\rho }_0\frac{\partial \overline{{\psi }_0}}{\partial \mathbf{E}}\right):\dot{\mathbf{E}}+\left({\tau }_{\mathrm{p}}-{\rho }_0\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}\right):{\stackrel{\Delta }{ϵ}}_{\mathrm{ep}}+{\tau }_{\mathrm{p}}:{\stackrel{\Delta }{ϵ}}_{\mathrm{p}}+$$

$$+{\rho }_0\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}:\left({ϵ}_{\mathrm{ep}}{\mathbf{l}}_{\mathrm{p}}+{\mathbf{l}}_{\mathrm{p}}^{\mathrm{T}}{ϵ}_{\mathrm{ep}}\right)+\sum _{i=1}^2\left({\tau }_{\mathrm{ov}}^i-{\rho }_0\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}\right):{\stackrel{\Delta }{ϵ}}_{\mathrm{ev}}^i+$$

$$+\sum _{i=1}^2\left({\tau }_{\mathrm{ov}}^i:{\stackrel{\Delta }{ϵ}}_{\mathrm{v}}^i\right)+\sum _{n=1}^2\left[{\rho }_0\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}:\left({ϵ}_{\mathrm{ev}}^i{\mathbf{l}}_{\mathrm{v}}^i+{\mathbf{l}}_{\mathrm{v}}^{i\mathrm{T}}{ϵ}_{\mathrm{ev}}^i\right)\right]\ge 0$$

$${\mathbf{S}}_{\mathrm{eq}}={\rho }_0\frac{\partial \overline{{\psi }_0}}{\partial \mathbf{E}}$$

$${\tau }_{\mathrm{p}}={\rho }_0\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}$$

$${\tau }_{\mathrm{ov}}^i={\rho }_0\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}$$

$$\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}:\left({ϵ}_{\mathrm{ep}}{\mathbf{l}}_{\mathrm{p}}+{\mathbf{l}}_{\mathrm{p}}^{\mathrm{T}}{ϵ}_{\mathrm{ep}}\right)=2{ϵ}_{\mathrm{ep}}\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}:{\stackrel{\Delta }{ϵ}}_{\mathrm{p}}$$

$$\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}:\left({ϵ}_{\mathrm{ev}}^i{\mathbf{l}}_{\mathrm{p}}^i+{\mathbf{l}}_{\mathrm{v}}^{i\mathrm{T}}{ϵ}_{\mathrm{ev}}^i\right)=2{ϵ}_{\mathrm{ev}}^i\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}:\stackrel{\Delta }{ϵ}{}_{\mathrm{v}}^i$$

$${\rho }_0\left(\mathbf{I}+2{ϵ}_{\mathrm{ep}}\right)\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}:{\stackrel{\Delta }{ϵ}}_{\mathrm{p}}+\sum _{i=1}^2{\rho }_0\left(\mathbf{I}+2{ϵ}_{\mathrm{ev}}^i\right)\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}:\stackrel{\Delta }{ϵ}{}_{\mathrm{v}}^i\ge 0$$

$${\stackrel{\Delta }{ϵ}}_{\mathrm{p}}:=c_{\mathrm{p}}{\rho }_0||\stackrel{\Delta }{ϵ}||\left(\mathbf{I}+2{ϵ}_{\mathrm{ep}}\right)\frac{\partial {\bar{\psi }}_{\mathrm{p}}}{\partial {ϵ}_{\mathrm{ep}}}$$

$$\stackrel{\Delta }{ϵ}{}_{\mathrm{v}}^i:=\frac{{\rho }_0}{{\eta }_{\mathrm{v}}^i}\left(\mathbf{I}+2{ϵ}_{\mathrm{ev}}^i\right)\frac{\partial {\bar{\psi }}_{\mathrm{v}}^i}{\partial {ϵ}_{\mathrm{ev}}^i}$$

$${\psi }_0=\frac{C_1}{\alpha }\left[e^{\alpha \left(I_1-\ln I_3-3\right)}\right]+D_2{\left(\ln I_3\right)}^2$$

$${\psi }_{\mathrm{p}}=\frac{C_{1p}}{{\alpha }_{\mathrm{p}}}\left[e^{{\alpha }_{\mathrm{p}}\left(I_1^{\mathrm{ep}}-\ln I_3^{\mathrm{ep}}-3\right)}\right]+D_{2p}{\left(\ln I_3^{\mathrm{ep}}\right)}^2$$

$${\psi }_{\mathrm{v}}^i=\frac{C_{1v}^i}{{\alpha }_{\mathrm{v}}^i}\left[e^{{\alpha }_{\mathrm{v}}^i\left(I_1^{\mathrm{ev}i}-\ln I_3^{\mathrm{ev}i}-3\right)}\right]+D_{2\mathrm{v}}^i{\left(\ln I_3^{\mathrm{ev}i}\right)}^2$$

$$I_1(\mathbf{C})=\mathrm{tr}(\mathbf{C})=\mathbf{C}:\mathbf{I};I_3(\mathbf{C})=\det \mathbf{C}=J^2$$

Full text

1 Introduction

Peat is considered a challenging natural material in geotechnical engineering practice. Its high water content, organic and fibrous nature cause difficulties in sampling and laboratory testing of peat materials. The large deformations [1, 2, 3], structural anisotropy [4], time-dependent behaviour [5] and great spatial variability in properties and fibrosity [6] observed in peat upon mechanical loading are indicative of the geomechanical complexities inherent in both experiments and numerical modelling. This study aims to find a rational ingress to the complex mechanical phenomena by setting up clear and plausible theories of the peat behaviour and propose a constitutive relationship for peat within the modern mechanics context adhering to principles established in the theory of materials.

The majority of experimental research on peat was focused on determining the undrained shear strength and the (long-term) time-dependent behaviour associated with drainage, i.e. consolidation, in one-dimensional (1D) settings. The undrained shear strength was analysed with Mohr-Coulomb friction-cohesion failure criteria. As peat structure and morphology are different from those of granular mineral soils, including the presence of fibres, the physical meanings of friction and cohesion are equivocal [7]. Corresponding to consolidation tests performed on peat, several models were developed or adopted from clay modelling, such as rheological models [8, 9, 10, 11], the $C_{\alpha}/C_{c}$ concept [12, 13, 14], isotache models [15, 16, 17], and elastic visco-plastic (EVP) models [18, 19, 20, 21, 22].

Early rheological models have reproduced some observed features which were proposed based on postulated peat consolidation mechanisms by considering the material’s micro-structure, such as the two-level structural model for fibrous peat by Berry et al. [9] and the rheological model proposed by Barden [8] which treated the volumetric and deviatoric stress response as a single process by employing a generalization in terms of the effective stress ratio. Due to the nonlinearity and mathematical complexity of the rheological models, curve-fitting techniques became a mainstream topic after the 1970s where two popular approaches prevailed. The $C_{\alpha}/C_{c}$ concept described a unique relationship between the index of secondary compression $C_{\alpha}$ and the compression index $C_{c}$ and led to a controversy on the uniqueness of the end of primary (EOP) consolidation as well as the so-called hypotheses A and B [23, 24, 13, 25, 26, 27, 28]. Besides, the 1D consolidation/compression models developed based on either time-line theory [29] or isotache theory [30] were adopted to model peats.

EVP models for peat simulations were adopted from the constitutive models developed for clay [18, 20, 21] based on the time-line concept and thus routinely adopted critical-state elastoplastic models for 1D small-strain conditions. An equivalent time was defined as the time needed to creep from a reference time line to the current state [18] and a unique limit time line existed in strain–effective stress space at infinite equivalent time [20]. By defining a reference time line for the current strain and effective stress, negative equivalent time could be obtained as a function of the strain–effective stress states, illustrating how heuristic many of these modelling approaches remained. The models developed based on the time-line concept [29, 31, 32, 20] directly took time as a variable in describing the rate-dependence. Such an approach causes difficulties in relating the time elapsed since loading to the intrinsically time-dependent or rather rate-dependent behaviour of a soil [15]. In most cases, using time as an explicit variable driving material behaviour is not advisable in constitutive theory of materials, cf. Haupt [33] and others. O’Loughlin [22] found that the EVP model overpredicted fibrous peat settlement.

The isotache model [15] proposed a unique relationship between effective stress, strain and strain rate by defining a linear relationship between Hencky strain and the natural logarithm of intrinsic time as well as a linear relationship between Hencky strain and the natural logarithm of effective stress. The Hencky strain measure was chosen also to account for the large strains experienced by peat. Time as a variable was replaced by strain rate in the isotache model, however, an intrinsic time as a material parameter needed to be defined. This type of models was constrained by its curve-fitting origin from tests at a constant level of stress in 1D conditions [34]. Albeit successful in modelling peat oedometer tests, the isotache model has difficulties in simulating additional loading conditions other than the original loading-unloading fitted curves, such as unloading relaxation [35].

General constitutive laws describing the stress–strain–time relationships under any possible loading conditions have been derived from mechanisms which could plausibly be inferred to cause the observed constitutive behaviour [34]. To account for the rate-dependence, the overstress theory [36, 37] and the nonstationary flow surface theory [38] are two perspectives which split rate-dependent elastoplasticity into rate-independent elastic and rate-dependent plastic parts. Neither of the theories are able to describe the relaxation process initiated within the yield surface due to the inherent combination of the viscous with plastic effects. In recent years, several elastoviscoplastic models were proposed for peat, such as the 1D small strain elastoviscoplastic model by Madachi and Gajo [39], the 3D small strain hyperplastic model by Boumezerane et al. [40] and the small strain kinematic bubble model by Boumezerane [41]. The models for peat reported in the literature suffered from a lack of experimental evidence—particularly regarding plastic yield criteria for peat—and consequently inherited certain concepts from the conventional small-strain approaches used in clay modelling.

This paper adopts work performed for hyperviscoelastic modelling of biological tissues [42, 43] and extends it by a plastic component, within a thermodynamically consistent framework for simulating finite-strain rate-dependent peat behaviour under deviatoric conditions. The rationale for starting from deviatoric conditions will be explained in Section 2. Deviating from conventional approaches to modelling plastic deformations of peat, a plastic flow rule without a yield criterion was adopted [44], see Section 3. A numerical implementation of the devised model was realized relying on features like the Kelvin mapping [45], nested Newton-Raphson iterations as well as implicit time integration and is described in Section 4. The viscoelastic parameters were calibrated from undrained triaxial relaxation tests and the elatoplastic parameters were calibrated from undrained loading-unloading test strain irrecoverability, see Section 6. An initial validation of the constitutive model was realized in undrained triaxial testing of an undisturbed peat at six strain rates with loading/unloading relaxations and is presented in Section 7. The article closes with a discussion and perspectives on future work in Sections 8 and 9.

2 General Methodology

This study adopts salient features of the material theory outlined by Haupt [33]. The total stress (${\mathbf{S}}$) will be additively decomposed into a rate-independent equilibrium (elastoplastic) part (${\mathbf{S}}_{\mathrm{eq}}$) and a rate-dependent (viscous) overstress (${\mathbf{S}}_{\mathrm{ov}}$) in Eq. (1). The term equilibrium refers to the thermodynamic state of the rate-dependent material in which all dissipative processes have ceased provided constant external conditions. The equilibrium stress is assumed to be a rate-independent functional of the process history, whereas the overstress is a rate-dependent function of the process history.

[TABLE]

The overstress ${\mathbf{S}}_{\mathrm{ov}}$ in Eq. (1) then has to disappear asymptotically for sufficiently slow processes or relax to zero during the course of any static continuation [33]. The equilibrium state test, therefore, can be achieved by termination relaxations or defined by tests carried out at sufficiently slow rates. This has been experimentally illustrated in undrained triaxial relaxation tests on clay [46]. Additionally, Graham et al. [26] found that the strain-rate influence on the disturbed shear strength of clays appeared to be independent of soil plasticity or stress history during laboratory consolidation, and the gradient of the shear stress–$\log$ (strain rate) relationship decreased markedly at low strain rates with a threshold strain rate of 0.2 %/hour for Haney clay [47] and about 0.05 %/hour for Drammen clay [48]. A similar concept was defined as static yield surface by Sheahan [49], where a minimum undrained strength was eventually reached regardless of further strain-rate reductions, i.e. an equilibrium state was reached after the termination of creep tests.

Summarising the observed geomechanical behaviour reported in the literature, undrained peat can be categorised as a rate-dependent material with thermodynamic equilibrium hysteresis, corresponding to a viscoplastic constitutive model. A hyperviscoelastic model was adopted motivated by the success of similar material models in the description of elastomeric [50, 44] and biological materials [42, 43]. Particularly the latter share some fundamental defining features with peat, such as a highly hydrated solid matrix being composed of fibre-reinforced organic matter. The hyperviscoelastic model was extended by incorporating plasticity to better account for the observed undrained peat behaviour. The proposed model was calibrated and validated with undrained triaxial tests of an undisturbed Sphagnum peat.

Quantifying and modelling the rate-dependence of fluid-filled porous media in general and peat in particular is complicated by the superposition of several non-linear effects. For example, the apparent viscoelasticity of such materials is due to both flow-dependent features caused by the momentum interaction between the solid matrix and the mobile pore fluid, and flow in-dependent phenomena associated with dissipative mechanisms in the solid matrix itself. An independent and physically based quantification of both contributions is practically only possible if first the flow-dependent features are removed from the picture by performing isochoric tests (fluid flow is primarily caused by volume changes of the pore space in a hydro-mechanical setting). Based on the resulting quantification of solid matrix rate dependence, non-isochoric tests such as oedometer tests can then be performed and modelled with a biphasic theory based on the effective stress concept to account for flow-dependent phenomena.

Hence, unconsolidated undrained triaxial tests were utilized here to calibrate the initial constitutive model on the merits of (1) the controllable boundary conditions; (2) avoiding the heterogeneity induced by consolidation [51] and the Mandel-Cryer effect [52]; and (3) avoiding rate-dependence of undrained peat resulting from the relative motion of pore water through the macro- and micro-porous structure. In a second step not presented here, the total stress constitutive model can be straightforwardly extended to an effective stress model as indicated in the previous paragraph. The laboratory testing plan was set up to allow the definition of a suitable equilibrium state for the peat material in undrained triaxial tests, to calibrate the rate-dependent and rate-independent components of the proposed model, and to provide additional experimental data for model validation.

Note, that by way of convention, the sign conventions of mechanics (stresses and strains positive in tension, pressures positive in compression) are used throughout the modelling section. In the experimental section, as well as in all graphs presenting modelling results, stresses and strains are plotted with positive values in compression (sign convention of soil mechanics).

3 Constitutive Model

A macroscopic phenomenological model is proposed which can suitably describe the observed rate-dependent behaviour of peat materials at an engineering scale. As such, the model does not account for microstructural phenomena explicitly. As the current focus is on the rate-dependence and plasticity of the porous structure itself (cf. Section 2), the pore pressure was not accounted for by an explicit hydro-mechanical coupling but was lumped into the material response by means of a low compressibility resulting in quasi-isochoric behaviour. Furthermore, the constitutive model was derived firmly based on continuum mechanical and material theoretical principles [33] to ensure a priori the thermodynamical consistency of the model results under arbitrary three-dimensional finite deformations.

3.1 Rheological model

The phenomenological rate-dependent constitutive model consisted of several parts, i.e. nonlinear elasticity, plasticity and rate-dependence. A parallel arrangement of a spring element, a spring-friction element and two Maxwell elements (Fig. 1) representing the elasticity, elasto-plasticity and visco-elasticity on two time scales, respectively, was chosen. This choice has successfully captured the complex behaviour of elastomeric [50, 44] and biological materials [42, 43].

To facilitate an intuitive understanding of this structure, consider first a uniaxial setting under small-strain conditions. The total stress decomposes, according to Fig. 1, into elastic, elasto-plastic and visco-elastic components:

[TABLE]

For each layer of the parallel structure aside from the purely elastic spring, the total strain can be decomposed into elastic and inelastic parts:

[TABLE]

3.2 3D Formulation

The three-dimensional finite-strain formulation of the constitutive model can be arrived at by a generalization of the uniaxial and small strain formulation [53, 54]. In the sequel, material isotropy and isothermal conditions are assumed. An additive decomposition of the total Second Piola-Kirchhoff stress ${\mathbf{S}}$ analogous to Eq. (2) is introduced in accordance with the rheological model:

[TABLE]

For the finite-strain formulation, the additive decomposition of the total strain in Eq. (3) corresponds to a multiplicative decomposition [55] of the deformation gradient, ${\mathbf{F}}$, into elastic and inelastic parts:

[TABLE]

We further introduce the Green-Lagrange strain tensor, ${\mathbf{E}}$, defined as

[TABLE]

with the right Cauchy-Green deformation tensor

[TABLE]

Only in the intermediate configurations induced by the split in Eq. (5), the total strain decomposes additively into purely elastic parts of Green-Lagrangian type and purely inelastic parts of Almansi-Euler type, i.e.

[TABLE]

with

[TABLE]

Objective strain rates in the current and intermediate configurations were obtained by Oldroyd rates via Lie-type time derivatives:

[TABLE]

where the superimposed dot denotes the material time derivative, and the inelastic velocity gradients are given by

[TABLE]

In the intermediate configurations, the total rates $\overset{\Delta}{{\mathbf{\epsilon}}}$ of the elastoplastic and viscoelastic components decompose additively into elastic and inelastic parts according to

[TABLE]

where $\overset{\Delta}{{\mathbf{\epsilon}}}$ decomposes into $\overset{\Delta}{{\mathbf{\epsilon}}}{}^{i}_{\mathrm{v}}$ and $\overset{\Delta}{{\mathbf{\epsilon}}}{}^{i}_{\mathrm{ev}}$ analogously with the help of ${\mathbf{l}}^{i}_{\mathrm{v}}$.

The concept of dual variables [53] ensures that certain scalar physically relevant measures, such as the stress power, remain invariant under configurational transformations:

[TABLE]

where ${\mathbf{\tau}}_{\mathrm{p}}$ and ${\mathbf{\tau}}^{i}_{\mathrm{ov}}$ are the corresponding stresses in the intermediate configurations of the plastic and the $i^{\mathrm{th}}$ viscoelastic second Piola-Kirchhoff stresses, i.e.

[TABLE]

In the context of the rheological model shown in Fig. 1, the isothermal Helmholtz free energy is the sum of the strain energies of all spring elements. Thus, the specific Helmholtz free energy is assumed to decompose analogous to the stresses:

[TABLE]

Thermodynamic consistency is ensured here by deriving the constitutive relations from the entropy inequality according to the Coleman-Noll procedure [56]. Demanding non-negative entropy production rates for an isothermal process, the Clausius-Duhem inequality reduces to the Clausius-Planck relation given in a Lagrangian setting by

[TABLE]

Insertion of the stress decomposition, Eq. (2) along with the constitutive assumption (22), yields

[TABLE]

where use has been made of the concept of dual variables, cf. Haupt [44], Görke et al. [42] and others. Standard arguments lead to the stress relations:

[TABLE]

Assuming $\bar{\psi}_{\mathrm{p}}$ and $\bar{\psi}^{i}_{\mathrm{v}}$ are isotropic tensor functions of their arguments ${\mathbf{\epsilon}}_{\mathrm{ep}}$ and ${\mathbf{\epsilon}}^{i}_{\mathrm{ev}}$, respectively, i.e. $\psi({\mathbf{A}})=\psi({\mathbf{Q}}{\mathbf{A}}{\mathbf{Q}}^{\mathrm{T}})\ \ \forall{\mathbf{Q}}\in\mathrm{SO(3)}$, we find

[TABLE]

leaving the residual dissipation inequality

[TABLE]

to be satisfied. Following Haupt [44], non-negativity is ensured by specifying the plastic flow rule and the viscous strain rate in the corresponding intermediate configurations with non-negative parameters $c_{\mathrm{p}}$ and $\eta^{i}_{\mathrm{v}}$:

[TABLE]

where $\hat{{\mathbf{C}}}_{\mathrm{ep}}={\mathbf{I}}+2{\mathbf{\epsilon}}_{\mathrm{ep}}$ and $\hat{{\mathbf{C}}}^{i}_{\mathrm{ev}}={\mathbf{I}}+2{\mathbf{\epsilon}}^{i}_{\mathrm{ev}}$ are the elastic deformation tensor in the corresponding intermediate configurations of right Cauchy-Green type. For the generalisation of the proposed constitutive relationships, the plastic and viscous flow rules are introduced as by proportionality between the inelastic strain rates and corresponding stresses. The information on the yielding of peat materials in undrained conditions in the literature is insufficient to motivate the choice of a particular yield surface and parameterise it sufficiently. To avoid parameterisation of yield surfaces and plastic potentials based on insufficient experimental data, and by reasonably assuming that peat exhibits plastic effects very early on during its loading response in undrained conditions, the present plasticity model requires neither a yield function nor a plastic potential and can be parameterized by one material parameter only. For simplification, the same Helmholtz free energy density functions, $\psi=\rho_{0}\bar{\psi}$ for a modified Neo-Hookean model [42], were adopted for all the spring elements in the constitutive model:

[TABLE]

where the material parameters $C_{1}$,$C_{\mathrm{1p}}$ and $C^{i}_{\mathrm{1v}}$ are elastic moduli affecting both the distortional and the volumetric behaviour; $D_{2}$, $D_{\mathrm{2p}}$ and $D^{i}_{\mathrm{2v}}$ quantify bulk behaviour of the material; and $\alpha$, $\alpha_{\mathrm{p}}$ and $\alpha^{i}_{\mathrm{v}}$ control the nonlinearity of the constitutive relationship; $I_{1}$, $I_{3}$ and $I_{1}^{\mathrm{ep}}$, $I_{3}^{\mathrm{ep}}$ and $I_{1}^{\mathrm{ev}i}$, $I_{3}^{\mathrm{ev}i}$ are the first and third principal invariants of the right Cauchy-Green type tensors ${\mathbf{C}}$, $\hat{{\mathbf{C}}}_{\mathrm{ep}}$ and $\hat{{\mathbf{C}}}^{i}_{\mathrm{ev}}$, respectively. Expressing the invariants by quantities defined in the reference configuration yields [50, 42]:

[TABLE]

With the Helmholtz free energy functions (34) and (35), and by pulling back the inelastic flow rules (31) and (32), the flow rules in the reference configuration read:

[TABLE]

and the total second Piola-Kirchhoff stress becomes

[TABLE]

4 Numerical Implementation

To allow the simulation of laboratory tests under the assumption of homogenous stress and deformation states, the material model was implemented into a numerical algorithm solving for mechanical equilibrium [57]. The proposed constitutive relations contain rate formulations which were discretised in time using an implicit backward Euler scheme. Kelvin mapping was employed to represent tensor coordinates by vectors (for second-order tensors) and matrices (for fourth-order tensors) while maintaining tensor character (as opposed to the Voigt mapping) [45]. Original tensors are denoted by bold symbols (${\mathbf{E}}$, $\boldsymbol{\sigma}$) while corresponding Kelvin-mapped quantities are designated by sans-serif italic font ($\mathsfbfit{E}$, $\mathsfbfit{\sigma}$).

The total Lagrangian formulation of the weak form and the global Newton-Raphson iteration scheme for non-linear material models in a full finite element context are elaborated in Bucher et al. [58], Zienkiewicz et al. [59] and Nagel et al. [45, 60]. Here, due to the homogeneity assumption, a simplified version could be used.

4.1 Global Newton and Local Stress Iterations

Nonlinearities were resolved using two nested Newton iteration algorithms for the global and local levels, respectively. The global problem is given by balancing internal and external stresses $\mathsfbfit{S}_{\mathrm{int}}=\mathsfbfit{S}_{\mathrm{ext}}$, i.e. by decreasing the norm of the residual

[TABLE]

to below a tolerance: $\|\mathsfbfit{r}_{\mathrm{g}}\|<\epsilon_{\mathrm{tol}}$. Corresponding to the displacement in a finite element setting, the primary unknown of the global problem here is the Green-Lagrange strain tensor $\mathsfbfit{E}$. The resulting Newton-Raphson scheme to solve for $\mathsfbfit{r}_{\mathrm{g}}=\mathsfbfit{0}$ introduces the global Jacobian (tangent) matrix

[TABLE]

The evolution equations (39) and (40) in the hyperviscoplastic model have to be integrated at each load step to calculate the increments of the inelastic internal tensor variables ${\mathbf{C}}_{\mathrm{p}}$ and ${\mathbf{C}}_{\mathrm{v}}^{i}$. The differential equations are discretised in time using the implicit Euler backward scheme formulated in (44) considering a general ordinary differential equation $\dot{y}=f(y)$ for a time step in the interval $[t,t+\Delta t]$:

[TABLE]

The time discretisation of the differential and algebraic equations necessary to integrate the stress increment can be written in the form of Eq. (45) by defining a residual vector $\mathsfbfit{r}$ and the state vector $\mathsfbfit{z}=[\mathsfbfit{S}^{\mathrm{T}},\mathrm{internal\ variables}^{\mathrm{T}}]^{\mathrm{T}}$, where $\mathsfbfit{E}^{i}$ is the strain tensor at the $i^{\mathrm{th}}$ global iteration and is considered as fixed during the local iterations

[TABLE]

A Taylor series expansion of (45) considering only linear terms in the neighbourhood of a given iterative solution of $\mathsfbfit{r}_{j}$ and $\partial\mathsfbfit{r}/\partial\mathsfbfit{z}|_{j}$ yields:

[TABLE]

Once the iteration converged, the use of the total differential of $\mathsfbfit{r}$ directly provides the consistent tangent matrix for the global iteration:

[TABLE]

In the resulting linear system, the first entry of the solution $\mathrm{d}\mathsfbfit{z}/\mathrm{d}\mathsfbfit{E}^{t+\Delta t}$:

[TABLE]

is the material contribution $\mathrm{d}\mathsfbfit{S}/\mathrm{d}\mathsfbfit{E}$ to the tangent matrix $\mathsfbfit{C}^{i}$ required in Eq. (43). This approach of two nested Newton-Raphson iterations enables quadratic convergence of both the local and the global problem due to an algorithmically consistent linearisation [61].

For the proposed hyperviscoplastic model, the components of the state vector is $\mathsfbfit{z}=[\widetilde{\mathsfbfit{S}}^{\mathrm{T}},\mathsfbfit{C}_{\mathrm{p}}^{\mathrm{T}},\mathsfbfit{C}^{1}_{\mathrm{v}}{}^{\mathrm{T}},\mathsfbfit{C}^{2}_{\mathrm{v}}{}^{\mathrm{T}}]^{\mathrm{T}}$, where $\widetilde{\mathsfbfit{S}}=\mathsfbfit{S}/C_{1}$ is a dimensionless stress measure. The residual vector, Jacobian and the global tangent solution of the hyperplastic model is presented in Appendix A.

The numerical implementation of the hyperviscoplastic model was successfully verified in four steps against specific analytical solutions of the elastic, viscoelastic, and elastoplastic model components under suitable test conditions [62].

5 Undrained Triaxial Tests

Undrained triaxial tests of an undisturbed Sphagnum peat were carried out to calibrate and validate the proposed finite-strain constitutive model. The model was calibrated in two steps, viz. equilibrium state and rate-dependent components, respectively. Correspondingly, suitable undrained triaxial loading-unloading equilibrium test needed to be determined experimentally, and tests at different loading rates were carried out. The remainder of this section details how a suitable equilibrium state test for the peat material was defined and details the experimental plan for model calibration.

The undisturbed peat samples were collected from Clara bog, Ireland, from beneath the groundwater table (sampling method and detailed properties, cf. Zhang and O’Kelly [63]). Some selective properties of the tested peat material are presented in Table 1.

A standard triaxial test set-up that included two pressure-volume controllers (GDS Instruments Ltd.) was used for undrained loading-unloading and relaxation tests on vertical undisturbed peat specimens ($\emptyset=38$ mm, $H=76$ mm). Saturation was confirmed by a Skempton $B$ value in excess of 0.95 for all specimens. The effect of cell pressure was investigated by undrained loading-unloading tests up to 5 % axial strain with incremental cell pressures of 0 kPa, 20 kPa, 40 kPa and 120 kPa on an undisturbed peat specimen at an axial strain rate of 16.0 %/hour. The test results demonstrated that the cell pressure had a negligible effect on the strain recovery as well as on the differential stress–axial strain relationships of the loading-unloading tests in this strain region. Thus, a confining pressure of 0 kPa was taken for all subsequent undrained triaxial tests.

5.1 Undrained triaxial relaxation tests

The transient behaviour of a soil is believed to be caused by viscous effects in the soil structure, resulting in a rate-dependent material response with observable features such as creep and stress relaxation. Since the viscous effects observed in soil deformation result from different material effects caused by intrinsic properties or phase interactions, the term "rate-dependent" provides a more accurate description than "time-dependent" and allows a wider perspective for interpretation. The rate-dependence of the tested undisturbed peat was studied in three steps for calibrating the proposed model: (1) two stress relaxation tests were carried out at two strain rates with pore water pressure measurement under undrained conditions; (2) comparison of a slower loading-unloading test to a loading-unloading test with step relaxations at a higher strain rate; (3) definition of equilibrium state test.

A triaxial relaxation test was carried out by stopping the displacement loading at a specific axial strain for a specimen tested at a given strain rate. Comparing to creep tests, which are commonly used for investigating material rate-dependence, relaxation tests record the stress change over time at a constant applied strain. Fig. 2 presents two triaxial compression tests to an axial strain of 5 % at axial strain rates of 160.0 %/hour and 16.0 %/hour followed by a 10 hour and a 24 hour relaxation phase, respectively, with pore water pressure measurement and zero cell pressure. The pore pressure variations were very small during the stress relaxations. A linear relationship can be found between the differential stress and the logarithm of time after an initial stress drop. Note that the "differential stress" is occasionally referred to as "deviator stress" in the literature. The latter terminology will be avoided here as the term "deviator" will be reserved for the actual deviatoric part of a tensor according to its mathematical definition which does not coincide with the differential stress.

It is reasonable to postulate the existence of relaxation termination points for the undrained relaxation tests, i.e. an undrained equilibrium state, as the stress dissipation rate reduced exponentially with relaxation duration. An undrained triaxial loading-unloading test setup was used to compare the results of a quicker test with step relaxation periods to a slower but monotonic test under otherwise identical testing conditions. Specifically, the undrained loading-unloading triaxial test performed at 160.0 %/hour with step relaxations for 16 min was compared with an undrained loading-unloading triaxial test at 4.81 %/hour. Fig. 3 illustrates that the stress–strain curve obtained by the slow monotonic test was in close vicinity of the connected stress-relaxation end points of the fast step-relaxation (160 %hour, 16 min) test. It was concluded that stress relaxation points of a higher-rate test with relaxations of a sufficient period of time can be approximated by a lower-rate test carried out at otherwise identical conditions. This observation was used for defining the equilibrium state.

5.2 Equilibrium state

The equilibrium state can be obtained by connecting the stress relaxation termination points or approximated by tests carried out at a low strain rate. To determine the termination of stress relaxation, an undrained stress relaxation test at an axial strain rate of 16.0 %/hour was relaxed at strain levels of 2 %, 4 % and 6 % for 24 hours, respectively. The three stress relaxation phases are plotted against time on normal scale axis in Fig. 4. The temperature influence on the stress relaxations was eminent by daily fluctuations. The temperature impact on the stress measurement turned out to be quite significant with an average sensitivity of 2.11 kPa/K making it difficult to precisely determine the relaxation termination points based on the current experimental equipment.

Connecting the final stress relaxation points at each strain level and comparing the resulting curve to the stress–strain curve obtained by a monotonic undrained compression test at an axial strain rate of 0.16 %/hour, again yielded a close proximity between the two curves as observed in the previous paragraph for different experimental settings (Fig. 5). Based on the finding that a lower-rate test can be used to approximate the relaxation termination points of a higher-rate test, we defined the undrained triaxial test at an axial strain rate of 0.16 %/hour as the “equilibrium test” for the undisturbed peat under undrained conditions for our current purposes.

By way of additional validation of this choice, a stress relaxation test was carried out at an axial strain rate of 0.16 %/hour to an axial strain of 6 %. Fig. 6 presents the stress relaxation against logarithmic time. The average peak differential stress dissipation rate (slope gradient of the steepest stress drop in Fig. 6) was determined to be 0.0038 kPa/min. Comparing this with the temperature sensitivity of the differential stress measured earlier, the highest amount of stress relaxation during a nine hour period was of the same magnitude as the stress variation instigated by a temperature change of 1.0 K. Therefore, the definition of the equilibrium state of the undisturbed peat specimens in undrained triaxial tests at an axial strain rate of 0.16 %/hour can be reasonably justified by both the relaxation termination curve as well as the stress dissipation rate during relaxation. No data has been reported in the literature regarding the non-stationary yield surface or equilibrium state testing for peat in either drained or undrained conditions. However, Vaid et al. [47] reported a threshold strain rate of rate-independence for undrained Haney Clay of 0.2 %/hour, which is in close proximity to the value chosen here.

Fig. 7 presents the differential stress – axial strain behaviour of the undisturbed peat at 0.16 %/hour. For a total axial strain of 20 %, the recovered strain based on the shape of the unloading curve was about 11.5 %. Fig. 7 indicates strain irrecoverability in the defined equilibrium state test after a full strain cycle (equilibrium hysteresis) of the undisturbed peat material under undrained testing conditions. The difference between the differential stress values from Fig. 7 and Fig. 5 was attributed to specimen variations as well as the temperature influence. Negative differential stress was recorded at the end of the unloading phase potentially resulting from the expulsion of the triaxial cell water (including the changes of hydrostaic pressure on the specimen and buoyancy on the load measurement) as well as inaccuracies due to temperature fluctuations.

Finally note, that the testing conditions characterizing the equilibrium state may be redefined, given more accurate experimental data. While this may improve the model fits quantitatively, no fundamental differences to the subsequent results are anticipated.

5.3 Undrained loading-unloading-reloading tests

In addition to the calibration tests described above, undrained loading-unloading-reloading triaxial tests at zero cell pressure were carried out at axial strain rates of 160 %/hr with step relaxations, and at 16.0 %/hr, 1.6 %/hr, 4.81 %/hr and 0.16 %/hr for model validation. The test results will be presented together with model validation in Section 7.

6 Parameter Fitting

The hyperviscoplastic model proposed for undrained peat behaviour has 15 material parameters in total, viz. 8 hyperelastic material parameters ($C_{1}$, $D_{2}$ for the elastic response; $C_{\mathrm{1p}}$, $D_{\mathrm{2p}}$ for the elasto-plastic response; $C_{\mathrm{1v}}^{i}$, $D_{\mathrm{2v}}^{i}$ where $i=1,\,2$ for the visco-elastic response); 4 parameters ($\alpha$, $\alpha_{\mathrm{p}}$, $\alpha_{\mathrm{v}}^{i}$) controlling the non-linearity of the constitutive relationship; one plastic parameter $c_{\mathrm{p}}$; and 2 viscosity parameters $\eta_{\mathrm{v}}^{i}$. For the undrained test conditions, the assumptions of incompressible and immiscible constituents of peat leads to a constant third principal invariant of the Right Cauchy Green tensor $I_{3}\approx 1.0$. In the numerical implementation of the constitutive model, material deformation is the solution of the global Newton-Raphson iterations. The strain-controlled load and cell pressure of 0 kPa represent the boundary conditions in the experimental triaxial tests performed. Near incompressibility is achieved numerically by setting the bulk modulus-like parameters $D_{2}$ to a very high value, thus penalizing volumetric deformation. The $D_{2}$ parameters were thus not subjected to fitting.

6.1 Experimental Data

To fit the numerical simulation to the experimental data, both need to be expressed in the same stress and strain spaces. For the experimental data, the stress measure was Cauchy stress ${\mathbf{\sigma}}$, i.e. force per current area, and the recorded strain was engineering strain $\varepsilon_{\mathrm{eng}}=\Delta h/h_{0}$; whereas the numerical data are expressed in terms of second Piola-Kirchhoff stress ${\mathbf{S}}$ and Green-Lagrange strain ${\mathbf{E}}$, respectively. The experimentally measured Cauchy stress in the specimen’s axial direction equals the differential stress ($q=\sigma_{11}-\sigma_{33}$), i.e. $\sigma_{11}=q$, given a cell pressure of $\sigma_{33}=0\,\mathrm{kPa}$. The deformation gradient for the undrained triaxial compression can be obtained from the vertical compression $\lambda=1+\varepsilon_{\mathrm{eng}}$ as

[TABLE]

The Cauchy stress and the second Piola-Kirchhoff stress are related via

[TABLE]

where $J$ is the determinant of the deformation gradient tensor, with a value of unity for the undrained test condition, cf. Eq. (49). Following this post-processing of the experimental data, model parameters were fitted by comparing the numerical simulation results with the experimental data.

6.2 Equilibrium State Test

Material parameters for the rate-independent part of the model were first fitted to the equilibrium test data. Thus the hyperviscoelastic components of the proposed model were eliminated by setting the $i^{\mathrm{th}}$ viscoelastic parameters to vanishing values. To facilitate the fitting, the material parameter sensitivities of the constitutive relationships were studied.

The experimental equilibrium test was simulated by prescribing boundary conditions as 0 kPa cell pressure and an axial strain rate of 0.16 %/hour compressing the test specimen to an axial engineering strain of 20 %, followed by unloading at the same rate. When the stress of the axial direction reached zero during the unloading, the strain-controlled load case switched to stress-control maintaining the vertically stress-free state. Material parameter sensitivity was tested with the hyperelastoplastic parameters listed in Table 2 and the results are plotted against the corresponding experimental results in Fig. 8.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Assigned with a relatively large value (500 kPa), parameters $D_{2}$ and $D_{\mathrm{2p}}$ had no influence on the constitutive relationship and successfully maintained the isochoric nature of the tests. Accordingly, the third invariant of the right Cauchy-Green deformation tensor in the simulation remained at 1.0$\pm$0.005. The loading curve was a hybrid effect of the parameters $C_{1}$, $C_{\mathrm{1p}}$, $c_{\mathrm{p}}$, $\alpha$ and $\alpha_{\mathrm{p}}$. A higher $C_{1}$ indicated a larger loading stiffness and a smaller irrecoverable plastic strain. $C_{\mathrm{1p}}$ determined the elastoplastic stress, thus influencing the plastic strain via the plastic flow rule. A larger $c_{\mathrm{p}}$ transformed a greater amount of the elastoplastic stress into plastic strain. The shape parameters $\alpha$ and $\alpha_{\mathrm{p}}$ had the same effect on the elastic spring and the elastoplastic spring elements, where $\alpha$ had a larger impact on the nonlinearity of the constitutive relationship with the selected material parameters.

The objective of the material parameter fitting was to find the proper combination of these parameters to balance between the loading curve and the irrecoverable plastic strain during unloading. The curve fitting process was carried out manually. Based on the hyperelastoplastic model parameter sensitivity tests, a set of the elastoplastic parameters was chosen that satisfactorily fulfilled this requirement. This set is presented in Table 3 with the model fitting results shown in Fig. 9.

6.3 Rate-dependent Tests

Freezing the fitted material parameters for the hyperelastoplastic model, tests at higher strain rates are needed for fitting the rate-dependent parts of the model. Undrained triaxial compression tests at two strain rates with relaxation were used for the parameter fitting of the two rate-dependent hyperviscoelastic model components. As outlined earlier, the parameters $D_{\mathrm{2v}}^{1}$ and $D_{\mathrm{2v}}^{2}$ have no influence on the simulated constitutive relationships under the present considerations, and were thus assigned with same values as $D_{2}$ and $D_{\mathrm{2p}}$ in the hyperelastoplastic component. Analogous to the hyperelastoplastic model, $\alpha_{\mathrm{v}}^{1}$ and $\alpha_{\mathrm{v}}^{2}$ are assigned to be zero, simplifying the Helmholtz free energy function (35) to a Neo-Hookean model. The model parameters $C_{\mathrm{1v}}$ and $\eta_{\mathrm{v}}$ of each viscoelastic component obtained by fitting the undrained triaxial compression tests with relaxation at the compression strain rates of 1.6 %/hour and 16.0 %/hour are given in Table 4. The results are presented in Figs. 10 and 11, respectively. Experimental and numerical results show some deviations but agree at least qualitatively to a satisfactory level considering the inter-specimen variations for peat materials and the temperature effect experienced during the experimental tests.

7 Model Validation

The hyperviscoplastic model with material parameters from Tables 3 and 4 was compared against undrained triaxial tests performed on five undisturbed vertical peat specimens at axial engineering strain rates of 160.0, 16.0, 4.81, 1.60 and 0.16 %/hour, which had not been used previously for fitting the model parameters. The strain-controlled loads from the experiments were directly applied to the numerical simulations. To compare the proposed model performance with the experimental results in both compression and unloading stages, the load cases for the following validations were in an utterly strain-controlled manner, thus tensile stresses were created during unloading as no platen lift-off was allowed.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The experimental data presented along with the model validation demonstrated the challenges in testing peat, viz. the intrinsic specimen variability, the manifest sensitivity on load measurements and temperature fluctuations due to the relatively low stress level. The equilibrium irrecoverable strains in the loading-unloading cycles to an axial strain of 20 % range from 2.8 % to 8.5 %. The hyperelastoplastic model parameters were obtained by fitting to the equilibrium loading-unloading test with a relatively large irrecoverable strain (8.5 %), thus resulting in a significant discrepancy between the simulation and the experimental results in unloading for the other tests which showed a quantitatively quite different hysteresis behaviour. Taking into account this large inter-specimen variability and the temperature effect, the simulation results in Fig. 12 show that the proposed hyperviscoplastic model was able to at least qualitatively predict peat behaviour in undrained triaxial compression tests for five different strain rates with step relaxations to which it was not fitted. This was true even for strain-rate ranges outside of the fitting domain. To test the unloading performance of the proposed model, the model was allowed to generate tensile stress at the end of unloading as the strain-controlled load was continuously applied albeit the stress reached negative values, whereas in the laboratory experiments, the load cell detached from the test specimens. Some initial remediation of this effect will be indicated Section 8.

8 Discussion

The validation of the hyperviscoplastic model demonstrated the model’s principal capability to predict other tests in a semi-quantitative fashion, with a better performance during the loading stage of the undrained triaxial tests. The relatively larger discrepancy during unloading between the numerical simulation and the experimental results could be improved by correcting the load case and modifying the plastic flow rule.

8.1 Load case modification

The load case modification was realized by allowing a detachment of the loading platen from the material. During the separation period between the loading platen and the material, the material was allowed to creep under a vanishing applied stress. When the loading platen reached the position of the creeping material during reloading, the strain-controlled load was reactivated. The hyperviscoplastic model with the modified load case is re-validated against the equilibrium test in Fig. 12 (e). Fig. 13 presents the modified load case which better reflected the test loading conditions.

As the proposed model was capable of simulating tests carried out at different strain rates, the triaxial loading-unloading-reloading test shown in Fig. 12 (c) was used as a rate-dependent representative test with strain-cycles for further investigation of the behaviour under the modified load case. The simulation result is presented in Fig. 14 and the unloading tensile stress was banished by the modified load case. The reloading curve was shifted to a larger irrecoverable plastic strain from the former unloading.

8.2 Plastic flow rule modification

The plastic flow rule was modified by eliminating the plastic strain evolution during passive unloading. When the applied load switched from loading to unloading, the material was still within the compressive domain. The original plastic flow rule in Eq. (39) determines the plastic strain increment direction as co-aligned with the current stress in the plastic element. Therefore, a modification of the flow rule was performed by suppressing plastic flow during passive unloading. Plastic flow appeared again during active load reversal. The modification was realized by multiplying a Heaviside step function to the current plastic flow rule that activates flow only if strain evolution and the stress in the plastic element driving the flow are aligned in the sense of a norm. Specifically, the current plastic flow rule (31) was modified in the intermediate configuration to

[TABLE]

where $H(\bullet)$ is the Heaviside step function and $\mathrm{sign}(\bullet)$ is the signum function.

The previously determined material parameters were used for the comparison of the hyperviscoplastic model with the modified plastic flow rule to the equilibrium test and the test at 1.6 %/hour. Fig. 15 and Fig. 16 present the simulation results of the tests carried out at 0.16 %/hour and 1.6 %/hour with the modified and the original flow rules. The partial plastic strains on the reference configuration, obtained by ${\mathbf{E}}_{\mathrm{p}}=\frac{1}{2}\left({\mathbf{C}}_{\mathrm{p}}-{\mathbf{I}}\right)$, were compared. The plastic strain increase during the unloading compressive domain was eliminated and remained constant instead. The irrecoverable plastic strain of the constitutive relationship during unloading was lessened by the modified plastic flow rule bringing it closer to the experimentally observed values.

9 Conclusions and Future Work

Peat is a challenging natural material for both laboratory testing and numerical modelling. Among the outstanding features of this article in the context of the literature available on peat modelling, are the systematic derivation of a finite-strain constitutive model within a modern continuum mechanical thermodynamically consistent framework in combination with modelling-oriented laboratory experiments designed to extract information on flow-independent viscoelastic phenomena. A hyperviscoplastic model with its elastic, viscoelastic,and elastoplastic components connected in parallel was adopted for modelling peat. The model provided qualitatively good simulation results of the undrained compression of the tested peat and demonstrated its versatility in simulating different load cases (loading-unloading with step relaxations at different strain rates) as well as its potential for further developments.

A certain advantage of a constitutive set-up like the presented one is its modular structure, with equilibrium and non-equilibrium parts modelled separately according to the experimental information available. Thus, certain "layers" of the model can simply be switched on or off depending on the modelling purpose and the availability of experimental data. The model can be improved in several ways, which are currently under study.

Further experimental data with better control of laboratory conditions are required for peat, such as keeping temperature fluctuations within $\pm 0.05$ K. Additionally, a statistical investigation of the possible spread of sample-specific equilibrium states for this type of peat material should be proposed based on a relatively large quantity of tests. Specifically, the definition of the equilibrium hysteresis must be defined in a manner that is consistent with the hysteresis observed in higher-rate tests. The strain-energy functions and flow-rule definitions chosen here were selected due to their generic nature and their low number of parameters. This makes the model versatile but leaves room to devise energy functionals and flow rules tailored to the characteristics of peat. Such a step requires further experimental study but is expected to greatly improve the model’s interpretive and predictive capabilities. Nevertheless, considering the generic nature of the expressions used here, the model already captured a wide spectrum of the experimental observations.

Furthermore, the model will be applied to peat under drained conditions. For that case, the various $D_{2}$-parameters governing the volumetric behaviour will have to be fitted in addition to quantifying flow-dependent effects by accounting for fluid-solid interaction. The model will, for that purpose, be included into a framework for hydromechanically coupled analyses, e.g. Nagel and Kelly [43] and Görke et al. [42], to capture transient effects linked to fluid pressurisation such as flow-dependent viscoelasticity. For that purpose, the devised model will constitute the constitutive law for the effective stresses ${\mathbf{S}}_{\mathrm{S}}^{\mathrm{E}}$

[TABLE]

where the pore pressure $p$ is derived from a saturation-constrained mass balance equation of the porous medium extended by a flow-model, e.g. Darcy’s law.

These extensions will include an implementation into a finite element framework to allow simulations with spatially heterogeneous solution and parameter fields. At that stage, the model will also be linked to advanced parameter identification algorithms to obtain material parameters by fitting to a series of tests simultaneously in an automated manner.

Acknowledgements

The authors acknowledge Mr. Eoin Dunne and Dr. Kevin Ryan for their assistance in collecting the peat samples as well as setting up the laboratory testing. This research was funded by Postgraduate Embark Initiative (grant number RS/2011/271) from the Irish Research Council and a Postgraduate Studentship Award from Trinity College Dublin. The modelling work was performed during a research stay at the Helmholtz-Centre for Environmental Research - UFZ. We cordially thank Prof. Olaf Kolditz for hosting the first author during that period, and the School of Engineering, Trinity College Dublin, for providing a travel and subsistence allowance for the same.

Appendix A Appendix: Local matrix-vector system

Residual vector of the numerical implementation for the hyperviscoplastic model:

[TABLE]

The state vector of the numerical implementation for the hyperviscoplastic model is:

[TABLE]

To solve the state vector, the Jacobian is of size $4\times 4$. The components of the hyperviscoplastic Jacobian are:

[TABLE]

Components of the RHS matrix for global tangent solution:

[TABLE]

[TABLE]

For the numerical implementation of the hyperviscoplastic model with modified plastic flow rule, three terms need to be updated in the above formulations:

Residual vector:

[TABLE]

Jacobian:

[TABLE]

RHS matrix for global tangent solution:

[TABLE]