## Introduction

The 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 Name | Description |
---|---|

muA | Shear modulus of network A |

thetaHat | Temperature factor |

lambdaL | Locking stretch |

kappa | Bulk modulus |

tauHatA | Flow resistance of Network A |

a | Pressure dependence of flow |

mA | Stress exponential of Network A |

n | Temperature exponetial |

muBi | Initial shear modulus of Network B |

muBf | Final shear modulus of Network B |

beta | Evolution rate of muB |

tauHatB | Flow resistance of Network B |

mB | Stress exponential of Network B |

muC | Shear modulus of Network C |

q | Relative contribution of I2 |

alpha | Thermal expansion coefficient |

theta0 | Thermal 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.

Index | State Variable Name |
---|---|

1 | Simulation time |

2 | Viscoelastic strain magnitude |

3 | Chain strain |

4 | Failure flag |

5-13 | Viscoelastic deformation gradient |

14-22 | Plastic deformation gradient |

23 | Shear modulus of network B |

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

Index | State Variable Name |
---|---|

1 | Simulation time |

2 | Viscoelastic strain magnitude |

3 | Chain strain |

4 | Failure flag |

5-10 | Viscoplastic Finger deformation tensor of Network A |

11-16 | Viscoplastic Finger deformation tensor of Network B |

17 | Shear modulus of network B |

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