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(double *stress, double *dStressDDStrain, const double *dStrain, const double *timeOld, const double dT, double &pNewDT)

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:

• stress – [inout] Cauchy stress

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

• dStrain – [in] linearized strain increment

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual int getNumberOfRequiredStateVars()

Returns:

Number of state variables required by the material.

virtual void assignStateVars(double *stateVars_, int nStateVars)

Assign state variable array to material.

Parameters:

• stateVars – [inout] Pointer to state variable array.

• nStateVars – [in] Number of state variables.

virtual StateView getStateView(const std::string &stateName)

Access material state variables by name.

Parameters:

stateName – [in] Name of the requested state variable.

Returns:

A view into the state variable array.

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].

std::unique_ptr<LinearViscoelasticPowerLawStateVarManager> stateVarManager

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

class LinearViscoelasticPowerLawStateVarManager : public MarmotStateVarVectorManager

Inheritance diagram for Marmot::Materials::LinearViscoelasticPowerLaw::LinearViscoelasticPowerLawStateVarManager:

Collaboration diagram for Marmot::Materials::LinearViscoelasticPowerLaw::LinearViscoelasticPowerLawStateVarManager:

Public Functions

inline LinearViscoelasticPowerLawStateVarManager(double *theStateVarVector, int nKelvinUnits)

Public Members

KelvinChain::mapStateVarMatrix kelvinStateVars

Public Static Attributes

static const auto layout = makeLayout( {{ .name = "kelvinStateVars", .length = 0 },} )