Linear Viscoelastic Power Law model

Theory

The present model is an implementation of an linear viscoelastic power law based on an approximation using a Kelvin chain. The constitutive law is given in total form as

\[\sig = \Cel : \epsE = \Cel : \left( \eps - \epsVE \right)\]

relating the nominal stress tensor \(\sig\) to the elastic strain tensor \(\epsE\) in terms of the isotropic fourth order stiffness tensor \(\Cel\). Where \(\eps\) is the total strain and \(\epsVE\) is the viscoelastic strain.

The rate of the compliance function of the material is given as

\[\JRate(t) = n m \frac{t^{n-1}}{\tau^n},\]

where \(n\) is the power law exponent, \(m\) is the power law compliance parameter, \(t\) is the time, and \(\tau\) is a reference time. The viscoelastic compliance is approximated by a Kelvin Chain.

The evolution of the viscoelastic strain \(\epsVE\) is defined as

\[\epsVERate(t) = \int_0^t \JRate(t-t') \DelNu : d\sig(t'),\]

in terms of the unit compliance tensor \(\DelNu\). The unit compliance tensor is defined as

\[\DelNu = \CelInv(E=1, \nu),\]

with the Young’s modulus \(E\) and the Poisson’s ratio \(\nu\).

Here we approximate a target power-law creep compliance \(J(t)\) by a Kelvin chain (generalized Kelvin–Voigt / Prony series) so the viscoelastic strain rate becomes a finite sum of exponentials that’s efficient to evaluate. The continuous form

\[J(t)=J_0+\int_0^\infty L(\tau)\,\bigl(1-e^{-t/\tau}\bigr)\,d\tau,\]

is replaced by

\[J_N(t)=J_0+\sum_{i=1}^N J_i\bigl(1-e^{-t/\tau_i}\bigr), \qquad \dot J_N(t)=\sum_{i=1}^N \frac{J_i}{\tau_i}\,e^{-t/\tau_i},\]

where \(J_i\) and \(\tau_i\) are the partial compliances and retardation times. For a specific compliance target the \(\{J_i,\tau_i\}\) are chosen so that \(J_N(t)\) tracks \(J(t)\) over the time window of interest, while the rate kernel \(\dot J_N\) approximates

\[\dot J(t)=\int_0^\infty \frac{L(\tau)}{\tau}\,e^{-t/\tau}\,d\tau.\]

The Post–Widder inversion formula provides a way to recover an approximate retardation spectrum from a known compliance function \(J(t)\).

\[L(\tau) \approx \frac{(-1)^{k}}{(k)!} \left(\frac{k}{\tau}\right)^{k} J^{(k)}\!\left(\frac{\tau}{k+1}\right), \qquad k \to \infty.\]

The order \(k\) defines the order of approximation of the actual continous spectrum. If it tends to infinity the continous spectrum is recovered. The discrete spectrum of the Kelvin chain is here used to calculate the stiffnesses or compliances for chosen retardation times by fitting the approximation of the continous spectrum.

A detailed description of the employed approach can be found, e.g., in Bazant & Jirasek (2018) Creep and Hygrothermal Effects in Concrete Structures.

Implementation

class LinearViscoelasticPowerLaw : public MarmotMaterialHypoElastic

Inheritance diagram for Marmot::Materials::LinearViscoelasticPowerLaw:

Collaboration diagram for Marmot::Materials::LinearViscoelasticPowerLaw:

Implementation of an isotropic linear viscoelastic material using a Power Law compliance function and assuming constant Poisson’s ratio generalized for 3D stress states.

Public Functions

LinearViscoelasticPowerLaw(const double *materialProperties, int nMaterialProperties, int materialLabel)

virtual void computeStress(state3D &state, Marmot::Matrix6d &dStressDDStrain, const Marmot::Vector6d &dStrain, const timeInfo &timeInfo) const override

For a given linearized strain increment \(\Delta\boldsymbol{\varepsilon}\) at the old and the current time, compute the Cauchy stress and the algorithmic tangent \(\frac{\partial\boldsymbol{\sigma}^{(n+1)}}{\partial\boldsymbol{\varepsilon}^{(n+1)}}\).

Parameters:

• state – [inout] A state3D instance carrying stress, strain energy, and state variables

• dStress_dStrain – [inout] Algorithmic tangent representing the derivative of the Cauchy stress tensor with respect to the linearized strain

• dStrain – [in] linearized strain increment

• timeInfo – Structure carrying the current (pseudo-)time and the (pseudo-)time increment

virtual double getDensity(const double *stateVars) const override

Get the mass density of the material.

Parameters:

stateVars – Pointer to the state variable array

Returns:

Mass density of the material

inline MarmotMaterialHypoElastic(const double *matProperties_, int nMaterialProperties_, int materialNumber_)

Constructs the material with a given set of material properties and an identifier.

Parameters:

• matProperties_ – [in] Pointer to the array of material properties.

• nMaterialProperties_ – [in] Number of entries in matProperties_.

• materialNumber_ – [in] Integer identifying this material instance.

Private Members

const double &E

Young’s modulus.

E represents the Young’s modulus for isotropic linear elasticity. It is a reference variable to materialProperties[0].

const double &nu

Poisson’s ratio.

nu represents Poisson’s ratio for isotropic linear elasticity. It is a reference variable to materialProperties[1].

const double &m

power law compliance parameter

m represents the power law compliance parameter. It is a reference variable to materialProperties[2].

const double &n

power law exponent

n represents the power law exponent. It is a reference variable to materialProperties[3].

const size_t nKelvin

number of Kelvin units to approximate the viscoelastic compliance

nKelvin represents the number of Kelvin units to approximate the power law function. It is a reference variable to materialProperties[4].

const double &minTau

minimal retardation time used in the viscoelastic Kelvin chain

minTau represents the minimal retardation time used in the Kelvin chain. It is a reference variable to materialProperties[5].

const double &timeToDays

ratio of simulation time to days

timeToDays represents the ratio of simulation time to days. It is a reference variable to materialProperties[6].

KelvinChain::Properties elasticModuli

Young’s modulus of the nKelvin Kelvin units.

KelvinChain::Properties retardationTimes

retardation times of the nKelvin Kelvin units

double zerothKelvinChainCompliance

compliance of the zeroth Kelvin unit

Private Static Attributes

static int powerLawApproximationOrder = 2