# PolyUMod Three Network (TN) Model

## Introduction

The PolyUMod Three Network (TN) model is a material model specifically developed for thermoplastic materials. It has many features that are similar to the hybrid model, but is designed to be more numerically efficient. The TN model is also a specialization of the more general Parallel Network model. This model is part of the PolyUMod library, and has Material Model (MM) id=11.

The TN model uses the following material parameters:

Parameter NameDescription
muAShear modulus of network A
thetaHatTemperature factor
lambdaLLocking stretch
kappaBulk modulus
tauHatAFlow resistance of Network A
aPressure dependence of flow
mAStress exponential of Network A
nTemperature exponetial
muBiInitial shear modulus of Network B
muBfFinal shear modulus of Network B
betaEvolution rate of muB
tauHatBFlow resistance of Network B
mBStress exponential of Network B
muCShear modulus of Network C
qRelative contribution of I2
alphaThermal expansion coefficient
theta0Thermal expansion reference temperature

#### General Notes about the Material Parameters:

• The material model can be made temperature independent by setting $$\hat{\theta}=0$$ and $$n=0$$.
• If temperature dependence is activated then the temperature should be in Kelvin or Rankine.
• There is no need to search for $$\lambda_L$$ unless some of the experimental data includes large strains.
• Set $$a=0$$ if only uniaxial tension or uniaxial compression data is available.
• The exponentials $$m_A$$ and $$m_B$$ should be less than 20 to ensure proper convergence.
• The flow resistance $$\hat{\tau}_A$$ should be less than $$\hat{\tau}_B$$.
• The parameter $$q$$ should be 0 unless biaxial experimental data is available.
• Set $$\alpha=0$$ during the calibration to prevent thermal expansion.

Here is an exemplar set of stress-strain predictions from the TN model created by MCalibration.

## TN Model in MCalibration

The material model can be solved using the internal MCalibration solver, or using any of the supported FE solvers. It is recommended to use the MCalibration native solver when possible since it runs significantly faster. If kappa=0, the applied load is uniaxial, and MCalibration is set as the solver, then a fast incompressible solver will be used to calculate the stress-strain response.

## Model Theory

As specified by its name, the kinematics of the three-network model consists of three parts, or molecular networks, acting in parallel, see the rheological representation figure to the right.

The total deformation gradient $$\mathbf{F}^{\mathit{appl}}$$ contains both a thermal expansion part $$\mathbf{F}^{\mathit{th}} = \left[1 + \alpha (\theta – \theta_0) \right] \mathbf{I}$$, and a mechanical deformation part $$\mathbf{F}$$:
$\mathbf{F}^{\mathit{appl}} = \mathbf{F}\, \mathbf{F}^{\mathit{th}}.$
The deformation gradient acting on network A is multiplicatively decomposed into elastic and viscoplastic components:
$\mathbf{F} = \mathbf{F}_A^e \mathbf{F}_A^v.$
The Cauchy stress acting on network A is given by a temperature-dependent version of the eight-chain representation:
$\boldsymbol{\sigma}_A = \frac{\mu_A} {J_A^e \overline{\lambda_A^{e*}}} \left[ 1 + \frac{\theta – \theta_0}{\hat{\theta}} \right] \frac{\mathcal{L}^{-1}\! \left( \overline{\lambda_A^{e*}} / \lambda_L \right) } {\mathcal{L}^{-1}\! \left( 1 / \lambda_L \right) } dev \left[ \mathbf{b}_A^{e*} \right] + \kappa (J_A^e – 1) \mathbf{I},$ where $$J_A^{e}=\det[\mathbf{F}_A^e]$$, $$\mu_A$$ is the initial shear modulus, $$\lambda_L$$ is the chain locking stretch, $$\theta$$ is the current temperature, $$\theta_0$$ is a reference temperature, $$\hat{\theta}$$ is a material parameter specifying the temperature response of the stiffness, $$\mathbf{b}_A^{e*} = (J_A^e)^{-2/3} \mathbf{F}_A^e (\mathbf{F}_A^e)^\top$$ is the Cauchy-Green deformation tensor, $$\overline{\lambda_A^{e*}} = \left(tr[\mathbf{b}_A^{e*}] / 3\right)^{1/2}$$ is the effective chain stretch based on the eight-chain topology assumption, $$\mathcal{L}^{-1}\!(x)$$ is the inverse Langevin function, where$$\mathcal{L}(x) = \coth(x)-1/x$$, and $$\kappa$$ is the bulk modulus. By explicitly incorporating the temperature dependence of the shear modulus it is possible to capture the stiffness variation of the material over a wide range of temperatures.

The viscoelastic deformation gradient acting on network B is decomposed into elastic and viscoplastic parts:
$\mathbf{F} = \mathbf{F}_B^e \mathbf{F}_B^v.$ The Cauchy stress acting on network B is obtained from the same eight-chain network representation that was used for network A.
$\boldsymbol{\sigma}_B = \frac{\mu_B} {J_B^e \overline{\lambda_B^{e*}}} \left[ 1 + \frac{\theta – \theta_0}{\hat{\theta}} \right] \frac{\mathcal{L}^{-1}\! \left( \overline{\lambda_B^{e*}} / \lambda_L \right) } {\mathcal{L}^{-1}\! \left( 1 / \lambda_L \right) } dev \left[ \mathbf{b}_B^{e*} \right] + \kappa (J_B^e – 1) \mathbf{1},$ where $$J_B^{e}=\det[\mathbf{F}_B^e]$$, $$\mu_B$$ is the initial shear modulus, $$\mathbf{b}_B^{e*} = (J_B^e)^{-2/3} \mathbf{F}_B^e (\mathbf{F}_B^e)^\top$$ is the Cauchy-Green deformation tensor, and $$\overline{\lambda_B^{e*}} = \left(tr[\mathbf{b}_B^{e*}] / 3\right)^{1/2}$$ is the effective chain stretch based on the eight-chain topology assumption.

The effective shear modulus is taken to evolve with plastic strain from an initial value of$$\mu_{Bi}$$ according to:
$\dot{\mu}_B = -\beta \left[ \mu_B – \mu_{\mathit{Bf}} \right] \cdot \dot{\gamma}_A,$ where$$\dot{\gamma}_A$$ is the viscoplastic flow rate. This equation enables the model to better capture the distributed yielding that is observed in many thermoplastics.

The Cauchy stress acting on network C is given by the eight-chain model with first order $$I_2$$ dependence:
$\boldsymbol{\sigma}_C = \frac{1}{1+q} \left\{ \frac{\mu_C} {J \overline{\lambda^*}} \left[ 1 + \frac{\theta – \theta_0}{\hat{\theta}} \right] \frac{\mathcal{L}^{-1}\! \left( \frac{\overline{\lambda^*}} { \lambda_L} \right) } {\mathcal{L}^{-1}\! \left( \frac{1} {\lambda_L} \right) } dev \left[ \mathbf{b}^* \right] + \kappa (J – 1) \mathbf{I} + q \frac{\mu_c}{J} \left[ I_1^* \mathbf{b}^* – \frac{2I_2^*}{3} \mathbf{I} – (\mathbf{b}^*)^2 \right] \right\},$ where $$J=\det[\mathbf{F}]$$, $$\mu_C$$ is the initial shear modulus, $$\mathbf{b}^* = J^{-2/3} \mathbf{F} (\mathbf{F})^\top$$ is the Cauchy-Green deformation tensor, and $$\overline{\lambda^*} = \left(tr[\mathbf{b}^*] / 3\right)^{1/2}$$ is the effective chain stretch based on the eight-chain topology assumption [Arruda:1993].

Using this framework, the total Cauchy stress in the system is given by $$\boldsymbol{\sigma} = \boldsymbol{\sigma}_A + \boldsymbol{\sigma}_B + \boldsymbol{\sigma}_C$$.

The total velocity gradient of network A,  $$\mathbf{L} = \dot{\mathbf{F}} \mathbf{F}^{-1}$$, can be decomposed into elastic and viscous components: $$\mathbf{L} = \mathbf{L}_A^e + \mathbf{F}_A^e \mathbf{L}_A^v \mathbf{F}_A^{e-1} = \mathbf{L}_A^e + \tilde{\mathbf{L}}_A^v$$, where $$\mathbf{L}_A^v = \dot{\mathbf{F}}_A^v \mathbf{F}_A^{v-1} = \mathbf{D}_A^v + \mathbf{W}_A^v$$ and $$\tilde{\mathbf{L}}_A^v = \tilde{\mathbf{D}}_A^v + \tilde{\mathbf{W}}_A^v$$. The unloading process relating the deformed state with the intermediate state is not uniquely defined since an arbitrary rigid body rotation of the intermediate state still leaves the state stress free. The intermediate state can be made unique in different ways, one particularly convenient way that is used here is to prescribe $$\tilde{\mathbf{W}}_A^v = \mathbf{0}$$. This will, in general, result in elastic and inelastic deformation gradients both containing rotations. The rate of viscoplastic flow of network A is constitutively prescribed by $$\tilde{\mathbf{D}}_A^v = \dot{\gamma}_A \mathbf{N}_A$$. The tensor $$\mathbf{N}_A$$ specifies the direction of the driving deviatoric stress of the relaxed configuration convected to the current configuration, and the term $$\dot{\gamma}_A$$ specifies the effective deviatoric flow rate. Noting that $$\boldsymbol{\sigma}_A$$ is computed in the loaded configuration, the driving deviatoric stress on the relaxed configuration convected to the current configuration is given by $$\boldsymbol{\sigma}_A’ = dev[\boldsymbol{\sigma}_A]$$, and by defining an effective stress by the Frobenius norm $$\mathbf{\tau}_A = || \boldsymbol{\sigma}_A’ ||_F \equiv \left( tr[\boldsymbol{\sigma}_A’ \boldsymbol{\sigma}_A’] \right)^{1/2}$$, the direction of the driving deviatoric stress becomes $$\mathbf{N}_A = \boldsymbol{\sigma}_A’ / \tau_A$$.
The effective deviatoric flow rate is given by the reptation-inspired equation [Bergstrom:2000]:
$\dot{\gamma}_A = \dot{\gamma}_0 \cdot \left(\frac{\tau_A}{\hat{\tau}_A + a R(p_A)} \right)^{m_A} \cdot \left(\frac{\theta}{\theta_0} \right)^n,$ where $$\dot{\gamma}_0 \equiv 1$$/s is a constant introduced for dimensional consistency, $$p_A = – [(\boldsymbol{\sigma}_A)_{11} + (\boldsymbol{\sigma}_A)_{22} + (\boldsymbol{\sigma}_A)_{33}]/3$$ is the hydrostatic pressure, $$R(x) = (x + |x|) / 2$$ is the ramp function, and $$\hat{\tau}_A$$,$$a$$,$$m_A$$, and $$n$$ are specified material parameters. In this framework, the temperature dependence of the flow rate is taken to follow a power law form. In summary, the velocity gradient of the viscoelastic flow of network A can be written
$\dot{\mathbf{F}}_A^v = \dot{\gamma}_A \mathbf{F}_A^{e-1} \frac{dev[\boldsymbol{\sigma}_A]}{\tau_A} \mathbf{F}.$

The total velocity gradient of network B can be obtained similarly to network A. Specifically, $$\mathbf{L} = \dot{\mathbf{F}} \mathbf{F}^{-1}$$ can be decomposed into elastic and viscous components: $$\mathbf{L} = \mathbf{L}_B^e + \mathbf{F}_B^e \mathbf{L}_B^v \mathbf{F}_B^{e-1} = \mathbf{L}_B^e + \tilde{\mathbf{L}}_B^v$$, where $$\mathbf{L}_B^v = \dot{\mathbf{F}}_B^v \mathbf{F}_B^{v-1} = \mathbf{D}_B^v + \mathbf{W}_B^v$$ and $$\tilde{\mathbf{L}}_B^v = \tilde{\mathbf{D}}_B^v + \tilde{\mathbf{W}}_B^v$$.

The unloading process relating the deformed state with the intermediate state is not uniquely defined since an arbitrary rigid body rotation of the intermediate state still leaves the state stress free. The intermediate state can be made unique in different ways [Boyce:1989], one particularly convenient way that is used here is to prescribe$$\tilde{\mathbf{W}}_B^v = \mathbf{0}$$. This will, in general, result in elastic and inelastic deformation gradients both containing rotations. The rate of viscoplastic flow of network B is constitutively prescribed by $$\tilde{\mathbf{D}}_B^v = \dot{\gamma}_B \mathbf{N}_B$$. The tensor $$\mathbf{N}_B$$ specifies the direction of the driving deviatoric stress of the relaxed configuration convected to the current configuration, and the term $$\dot{\gamma}_B$$ specifies the effective deviatoric flow rate. Noting that $$\boldsymbol{\sigma}_B$$ is computed in the loaded configuration, the driving deviatoric stress on the relaxed configuration convected to the current configuration is given by $$\boldsymbol{\sigma}_B’ = dev[\boldsymbol{\sigma}_B]$$, and by defining an effective stress by the Frobenius norm $$\mathbf{\tau}_B = || \boldsymbol{\sigma}_B’ ||_F \equiv \left( tr[\boldsymbol{\sigma}_B’ \boldsymbol{\sigma}_B’] \right)^{1/2}$$, the direction of the driving deviatoric stress becomes$$\mathbf{N}_B = \boldsymbol{\sigma}_B’ / \tau_B$$. The effective deviatoric flow rate is given by the reptation-inspired equation [Bergstrom:2000]:
$\dot{\gamma}_B = \dot{\gamma}_0 \cdot \left(\frac{\tau_B}{\hat{\tau}_B + a R(p_B)} \right)^{m_B} \cdot \left(\frac{\theta}{\theta_0} \right)^n,$ where $$\dot{\gamma}_0 \equiv 1$$/s is a constant introduced for dimensional consistency, $$p_B = – [(\boldsymbol{\sigma}_B)_{11} + (\boldsymbol{\sigma}_B)_{22} + (\boldsymbol{\sigma}_B)_{33}]/3$$ is the hydrostatic pressure, and $$\hat{\tau}_B$$,$$a$$,$$m_B$$, and $$n$$ are specified material parameters. In this framework, the temperature dependence of the flow rate is taken to follow a power law form. In summary, the velocity gradient of the viscoelastic flow of network B can be written
$\dot{\mathbf{F}}_B^v = \dot{\gamma}_B \mathbf{F}_B^{e-1} \frac{dev[\boldsymbol{\sigma}_B]}{\tau_B} \mathbf{F}.$

## Notes

There are currently two implementations of the TN model. The default implementation is used by all solvers except Abaqus/Explicit. The implementation that is used by Abaqus/Explicit is a newer implementation that is more numerically efficient, but only supports one choice of ODE solver and uses a different set of state variables. By setting global material parameter 1 to -11 (instead of 11), Abaqus/Explicit will use the default implementation instead of the new implementation.

## State Variables

The state variables that are used by the TN model for all FE solvers except Abaqus/Explicit are summarized in the following table.

IndexState Variable Name
1Simulation time
2Viscoelastic strain magnitude
3Chain strain
4Failure flag
23Shear modulus of network B

The state variables that are used by Abaqus/Explicit are summarized in the following table:

IndexState Variable Name
1Simulation time
2Viscoelastic strain magnitude
3Chain strain
4Failure flag
5-10Viscoplastic Finger deformation tensor of Network A
11-16Viscoplastic Finger deformation tensor of Network B
17Shear modulus of network B

Note that the LS-DYNA Explicit implementation is using 44 state variables.

### Convert Prony Series to Frequency Data

Instructions for how to quickly convert a Prony series to storage and loss modulus as a function of frequency.

### Convergence Properties for Different MCalibration Optimization Methods

MCalibration supports many different optimization methods. This study examines the performance of the different methods.

### Material Model for Polyethylene Terephthalate (PET)

Calibrated non-linear viscoplastic material model for PET.

### 4 thoughts on “PolyUMod Three Network (TN) Model”

1. I want to apply this TN model for PTFE parts in the assembly.
Can you please suggest the experiments to conduct. Are there any guidelines?

1. Thank you, Jorgen
I understand cyclic test can quickly help instead of multiple tests, and also capture relaxation behavior.
Uniaxial Compression and Tension are sufficient to calibrate the best fit models.
I will use Cyclic test-OPTION 1 to carry out further tests.
Wondering if there is a thumb rule for selecting
1. strain amplitudes,
2. no. of strain steps/hysteresis loops
3. relaxation time (how many seconds)

2. maximprincen

I can’t really figure out what these state variables are and what they do. Could you please elaborate on this?