# Three Network Model

The Three Network (TN) model was developed for predicting the large-strain, viscoplastic reponse of generic thermoplastic materials (in a glassy state). The TN-model is an excellent generic model for predicting the response of many different classes of both amorphous and semicrystalline thermoplastics.

## Parameters for the TN Model

IndexSymbolNameUnitDescription
1$$\mu_A$$muAStressShear modulus of network A
2$$\hat{\theta}$$thetaHatTemperatureTemperature factor
3$$\lambda_L$$lambdaL-Locking stretch
4$$\kappa$$kappaStressBulk modulus
5$$\hat{\tau}_A$$tauHatAStressFlow resistance of network A
6$$a$$a-Pressure dependence of flow
7$$m_a$$mA-Stress exponent of network A
8$$n$$n-Temperature exponential
9$$\mu_{Bi}$$muBiStressInitial shear modulus of network B
10$$\mu_{Bf}$$muBfStressFinal shear modulus of network B
11$$\beta$$beta-Evolution rate of muB
12$$\hat{\tau}_B$$tauHatBStressFlow resistance of network B
13$$m_B$$mB-Stress exponent of network B
14$$\mu_C$$muCStressShear modulus of network C
15$$q$$q-Relative contribution of I2 of network C
16$$\alpha$$alpha1/TemperatureThermal expansion coefficient
17$$\theta_0$$theta0TemperatureThermal expansion reference temperature

# Notes

• The temperature unit should be either Kelvin or Rankine when using a temperature-dependent TN model.
• The temperature dependence can be turned off by setting $$\hat{\theta}=0$$, and n=0.
• Set $$\alpha=0$$ if thermal expansion should not be included in the model calculations.
• If q = 0 then the material model can be calibrated using only uniaxial data. If q > 0 then experimental data from multiple loading modes should be used for the calibration.
• It is recommended to set $$\hat{\tau}_A \le \hat{\tau}_B$$, and that $$\mu_{Bi} \ge \mu_{Bf}$$.
• There is no need to search for $$\lambda_L$$ when calibrating the TN model unless the experimental data set includes data to very large strains.
• Set $$a=0$$ if only uniaxial tension or compression data is available.
• The exponentials $$m_A$$ and $$m_B$$ should be less than 20 to ensure proper convergence.
• The PolyUMod library contains 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  optimized for explicit calculations, but only supports one choice of ODE solver and uses a different set of state variables. By setting the global parameter 1 to a value of -11 (instead of 11), Abaqus/Explicit will use the default implementation instead of the optimized implementation.

# Model Theory

The details of the TN-model are described in the Mechanics of Solid Polymers book by J. Bergstrom. This section only summarizes the key features of the model theory.

As specified by its name, the kinematics of the Three-Network (TN) model consists of three parts, or molecular networks, acting in parallel, see the rheological representation 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}^{\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) } \text{dev} \left[ \mathbf{b}_A^{e*} \right] + \kappa (J_A^e – 1) \mathbf{1},$$
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(\text{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. Note if $$\hat{\theta}=0$$ then temperature dependence of the stiffness is disabled.

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) } \text{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(\text{tr}[\mathbf{b}_B^{e*}] / 3\right)^{1/2}$$ is the effective chain stretch based on the eight-chain topology assumption. In this equation  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 defined below. This equation enables the model to better capture the distributed yielding that is observed in many thermoplastics.

Note, to make this feature work as intended the flow of network A should occur before the flow of network B. That is, $$\hat{\tau}_A$$ should be less than $$\hat{\tau}_B$$. Also, $$\mu_{Bi}$$ should be larger than $$\mu_{Bf}$$.

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 \lambda_L} \left[ 1 + \frac{\theta – \theta_0}{\hat{\theta}} \right] \frac{\mathcal{L}^{-1} \left( \frac{\lambda_L} { \lambda_L} \right) } {\mathcal{L}^{-1} \left( \frac{1} {\lambda_L} \right) } \text{dev} \left[ \mathbf{b}^* \right] + \kappa (J – 1) \mathbf{1} + 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 $$\lambda_L = \left(\text{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 [Boyce, 1989], 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’ = \text{dev}[\boldsymbol{\sigma}_A]$$, and by defining an effective stress by the Frobenius norm $$\mathbf{\tau}_A = || \boldsymbol{\sigma}_A’ ||_F \equiv \left( \text{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{\text{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’ = \text{dev}[\boldsymbol{\sigma}_B]$$, and by defining an effective stress by the Frobenius norm $$\mathbf{\tau}_B = || \boldsymbol{\sigma}_B’ ||_F \equiv \left( \text{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{\text{dev}[\boldsymbol{\sigma}_B]}{\tau_B} \mathbf{F}.$$

# State Variables

The TN model uses the following state variables for all FE solvers except Abaqus/Explicit:

• State variable 1: Simulation time
• State variable 2: Viscoplastic strain magnitude
• State variable 3: Chain strain
• State variable 4: Failure flag
• State variables 5-13: Deformation gradient $$\mathbf{F}_A^v$$
• State variables 14-22: Deformation gradient $$\mathbf{F}_B^v$$
• State variable 23: Shear modulus $$\mu_B$$

The TN model for Abaqus/Explicit uses the following state variables:

• State variable 1: Simulation time
• State variable 2: Viscoplastic strain magnitude
• State variable 3: Chain strain
• State variable 4: Failure flag
• State variables 5-10: Tensor $$\mathbf{c}_{Av}^{-1}$$
• State variables 11-16: Tensor $$\mathbf{c}_{Bv}^{-1}$$
• State variable 17: Shear modulus $$\mu_B$$

# Exemplar Model Predictions

The following figure shows exemplar stress-strain predictions from the TN-model is uniaxial tension at different strain rates.