A Multiscale Multisurface Constitutive Model for The Thermo-Plastic Behavior of Polyethylene

TL;DR

This paper develops a multiscale model that predicts the thermo-plastic behavior of polyethylene by integrating molecular dynamics simulations with continuum mechanics, enabling accurate yield predictions across different strain rates and temperatures.

Contribution

The paper introduces a novel hierarchical multiscale modeling approach that links molecular simulations to continuum predictions for polyethylene's yield behavior.

Findings

Multisurface yield functions for polyethylene identified.

High strain rate MD simulations scaled to quasi-static conditions.

Temperature and strain rate dependent yield strength predicted.

Abstract

We present a multiscale model bridging length and time scales from molecular to continuum levels with the objective of predicting the yield behavior of amorphous glassy polyethylene (PE). Constitutive parameters are obtained from molecular dynamics (MD) simulations, decreasing the requirement for ad-hoc experiments. Consequently, we achieve: (1) the identification of multisurface yield functions; (2) the high strain rate involved in MD simulations is upscaled to continuum via quasi-static simulations. Validation demonstrates that the entire multisurface yield functions can be scaled to quasi-static rates where the yield stresses are possibly predicted by a proposed scaling law; (3) a hierarchical multiscale model is constructed to predict temperature and strain rate dependent yield strength of the PE.

Figures

Tables

Table 1. Table 1: Functional form and parameters of the Dreiding force field

Interaction	Form	Parameters
Bond	$E_b=\frac{1}{2}{k}_b{\left(r-r_{eq}\right)}^2$	$k_b=350kcal/mol{\AA }^2,r_{eq}=1.53{\AA }^2$
Angle	$E_a=\frac{1}{2}{k}_{\theta }{\left(cos\left(\theta \right)-cos\left({\theta }_{eq}\right)\right)}^2$	$k_{\theta }=60kcal/mol/ra{d}^2,{\theta }_{eq}={109.5}^0$
Dihedral	$E_d=\frac{1}{2}{\sum }_{i=0}^3{d}_{i}co{s}^i\left(\varphi \right)$	$d_0=1.736,d_1=-4.490,d_2=0.776,d_3=6.990\left(kcal/mol\right)$
Non-bonded		$\xi =0.112kcal/mol,\delta =4.01\AA ,r_{cut}=10.5\AA$

Table 2. Table 2: Multisurface constitutive model algorithm overview for PE.

(1)	Compute trial stress,
	${𝝈}_{n+1}^{tr}={𝝈}_n+{𝑫}^e:\Delta 𝜺$
	representing the stress in terms of the von-Mises equivalent $q_{n+1}^{tr}$ and hydrostatic $p_{n+1}^{tr}$ stresses:
	$q_{n+1}^{tr}=\sqrt{\frac{3}{2}{𝒔}_{n+1}^{tr}:{𝒔}_{n+1}^{tr}}$, $p_{n+1}^{tr}=p_n+K\Delta {\epsilon }_v$ with ${𝒔}_{n+1}^{tr}={𝒔}_n+2G\Delta {𝜺}_d$
(2)	Consider the temperature dependent elastic and yield behavior:
	The temperature ($T$) dependence of the Young’s modulus is explained by Equation (3):
	$E\left(T\right)=E^{ref}\left(log{a}_T\right)$
	The yield stresses and hardening laws dependent on the temperature ($T$)
	is described by Equation (19):
	${\sigma }_n={\sigma }_n^{ref}+{\beta }_n\left(T-T^{ref}\right)$
(3)	Check yield criterion given by Equation (11):
	IF $f(p,q,{\epsilon }_e^p)\le 0$ THEN
	${𝝈}_{n+1}={𝝈}_{n+1}^{tr}$, $q_{n+1}=q_{n+1}^{tr}$, $p_{n+1}=p_{n+1}^{tr}$ and EXIT
	ELSE
	Perform return mapping algorithm to obtain plastic multiplier $\Delta \lambda$.
	ENDIF
(4)	Update stress tensor
	$q_{n+1}=\sqrt{\frac{3J_2^{tr}}{{\left(1+6G\Delta \lambda \right)}^2}}$, ${𝒔}_{n+1}=\frac{{𝒔}_{n+1}^{tr}}{1+6G\Delta \lambda }$, $p_{n+1}=\frac{p_{n+1}^{tr}}{1+2K\alpha \Delta \lambda }$
	${𝝈}_{n+1}={𝝈}_{n+1}^{tr}-6G\Delta \lambda {𝒔}_{n+1}-\frac{2}{3}K\alpha \Delta \lambda p_{n+1}𝑰$
(5)	EXIT

Table 3. Table 3: Validation of the predicted tensile and compressive yield stresses with experimental result at different strain rates.

Deformation	Strain rate	Quasi-static simulations	Experimental results
Tension	0.2 $mi{n}^{-1}$	25.92 (MPa)	25.0 (MPa) [23]
	8.0 $mi{n}^{-1}$	28.84 (MPa)	30.8 (MPa) [23]
Compression		- 37.61 (MPa)	-31.72 (MPa) [29]

Table 4. Table 4: Constitutive properties for the PE model obtained from MD simulations at the room temperature.

$E$	$\nu$	${\nu }_p$	${\sigma }_t^{{\text{static}}^{ref}}$	${\sigma }_c^{{\text{static}}^{ref}}$	${\sigma }_s^{{\text{static}}^{ref}}$	${\beta }_T$	${\beta }_C$
1.32 GPa	0.32	0.32	25.92 MPa	-37.61 MPa	16.34 MPa	0.55	0.93
${\beta }_S$	$C_1$	$C_2$	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$
0.58	0.44	82.5	-0.66	-0.35	0.01	6.04

Equations

$${\mathbf{\epsilon }}_e=\sqrt{\frac{1}{2}}\sqrt{({\epsilon }_{11}-{\epsilon }_{22}{)}^2+({\epsilon }_{22}-{\epsilon }_{33}{)}^2+({\epsilon }_{11}-{\epsilon }_{33}{)}^2+\frac{4}{3}{\gamma }_{12}^2+\frac{4}{3}{\gamma }_{23}^2+\frac{4}{3}{\gamma }_{13}^2},$$

$${\dot{\mathbf{\epsilon }}}_e=\sqrt{\frac{1}{2}}\sqrt{({\dot{\epsilon }}_{11}-{\dot{\epsilon }}_{22}{)}^2+({\dot{\epsilon }}_{22}-{\dot{\epsilon }}_{33}{)}^2+({\dot{\epsilon }}_{11}-{\dot{\epsilon }}_{33}{)}^2+\frac{4}{3}{\dot{\gamma }}_{12}^2+\frac{4}{3}{\dot{\gamma }}_{23}^2+\frac{4}{3}{\dot{\gamma }}_{13}^2},$$

$$log{a}_T=\frac{-C_1(T-T^{ref}-140)}{C_2+(T-T^{ref}-140)}$$

$$\sigma =\frac{\Delta H}{V^{\ast }}+\frac{RT}{V^{\ast }}\ln \frac{2\dot{\gamma }}{{\dot{\gamma }}_0}$$

$${{\sigma }_t^{\text{static}}|}_{{}_{T^{ref}}}=\frac{{{\sigma }_t|}_{{}_{T^{ref}}}}{\gamma }\approx 25.92\text{MPa}\Rightarrow \gamma \approx 0.23$$

$${{\sigma }_c^{\text{static}}|}_{{}_{T^{ref}}}=\frac{{{\sigma }_c|}_{{}_{T^{ref}}}}{\gamma }\approx -37.6\text{MPa};{{\sigma }_s^{\text{static}}|}_{{}_{T^{ref}}}=\frac{{{\sigma }_s|}_{{}_{T^{ref}}}}{\gamma }a\approx 16.34\text{MPa}$$

$$p(\mathbf{\theta }|\mathbf{z})=\frac{p(\mathbf{z}|\mathbf{\theta })p(\mathbf{\theta })}{p(\mathbf{z})},$$

$$\underset{posterior}{\underbrace{p(\mathbf{\theta }|\mathbf{z})}}\propto \underset{likelihood}{\underbrace{p(\mathbf{z}|\mathbf{\theta })}}\underset{prior}{\underbrace{p(\mathbf{\theta })}}.$$

$${\mathbf{\theta }}_{MAP}=\underset{\mathbf{\theta }}{\mathrm{argmax}}p(\mathbf{z}|\mathbf{\theta })p({\mathbf{\theta }}_i)$$

$${\sigma }_t={\theta }_1{e}^{{\theta }_2\dot{\epsilon }}+{\theta }_3{e}^{{\theta }_4\dot{\epsilon }}$$

$$f(p,q,{\epsilon }_e^p)=q-\beta ({\epsilon }_e^p)p-c({\epsilon }_e^p)$$

$$\beta ({\epsilon }_e^p)$$

$$g=q^2+\alpha p^2$$

$$\alpha =\frac{9}{2}\frac{1-2{\nu }_p}{1+{\nu }_p}$$

$$\Delta {\epsilon }_{22}^p=\Delta {\epsilon }_{33}^p=-{\nu }_p\Delta {\epsilon }_{11}^p$$

$$\Delta {\mathbf{\epsilon }}^p=\Delta \lambda \frac{\partial g}{\partial \mathbf{\sigma }}$$

$${\epsilon }_e^p=\sqrt{k{\epsilon }^p:{\epsilon }^p}$$

$${\sigma }_t={\sigma }_t\left({\epsilon }_e^p\right),{\sigma }_c={\sigma }_c\left({\epsilon }_e^p\right),{\sigma }_s={\sigma }_s\left({\epsilon }_e^p\right),{\sigma }_{bt}={\sigma }_{bt}\left({\epsilon }_e^p\right),{\sigma }_{bc}={\sigma }_{bc}\left({\epsilon }_e^p\right)$$

$${\sigma }_n={\sigma }_n^{ref}+{\beta }_n\left(T-T^{ref}\right)$$

Full text

1 Introduction

Polymers are extensively used in industrial applications, particularly in the aerospace and automotive industry due to their physical and mechanical properties. Polymeric composites are highly complex and their mechanical properties depend on many variables such as temperature, strain rates, etc. and the physics of plastic deformations in amorphous polymers has not been well understood. Attempts were made to understand, especially, temperature and strain rate dependent yielding in polymers, see [1, 2] for a concise review. Nevertheless, these models tend to rely on experiments whose constitutive parameters are not physically motivated and can only be used to predict behavior of a specific material. In engineering practice, a visco-plastic model based on the pressure-modified von- Mises criterion is commonly used for thermoplastic polymers [3, 4]. However, these polymers behave differently under tensile, compressive and shear deformations. Hence, the von-Mises yield criterion is no longer appropriate. A number of theoretical studies have been done to find proper yield surface for the prediction of thermoplastics [4, 5]. Multisurface yield functions seem a suitable candidate to describe the yield behavior for a wide range of polymers [4]. However, calibration (fitting) procedures are not always possible for multiaxial loading conditions especially, due to loss of data.

To accurately predict macroscopic properties the molecular feature associated with the plastic mechanism must be understood [6]. Molecular theories of plastic behavior in amorphous polymers were reviewed by Stachurski [7]. However, the behavior at nano length scales was not explained within scope of these theories.

Along with the development of accurate inter-atomic potentials using quantum mechanics, molecular dynamics (MD) simulations are a powerful tool in visualizing molecular mechanisms of yielding [8]. MD simulations offer a promising way to develop new theories and models for glassy polymers as they can reduce the need for ad-hoc experiments.

Fully atomistic models based on force fields and chemical structure of materials allow us to physically interpret their complex physical phenomena [9], the length and time scales, nevertheless, limit the mechanism associated with viscoelastic/plastic behavior of the material, since interactions between single atoms are explicitly considered. Coarse graining methods, such as united atom (UA) models, can increase the length scales though the time scale has been still limited [6]. Rottler et al. [3] studied shear yielding of glassy polymers under multiaxial loading conditions using MD simulations and relate them to the pressure-modified von-Mises criterion, but just at nanoscale as the simulation time is prohibitive.

However, macroscopic continuum mechanics models can be employed to study large domains and realistic time-scales. Therefore, a multiscale model passing the nanoscale descriptions to the continuum is very important. In other words, macroscopic constitutive parameters describing the evolution of macroscopic properties can be obtained from MD simulations [10].

Due to high strain rates involved in MD simulations, which are not experimentally encountered, an appropriate scaling law for the yield surface [10] is essential to reconcile the different strain rates of MD simulations and experiments. In this article, the quasi-static simulations are employed to extract yield stresses at quasi-static strain rates from MD simulations. Furthermore, Bayes’ theorem is used to construct an upscaling technique. In particular, Bayesian approach considering prior information of the parameters (i.e. strain rate) upon which the posterior distribution is updated given a set of observations, leading to an identification of constitutive parameters.

The article begins with the nanoscale model of the PE. Temperature dependence of the elastic and yield behavior is subsequently accounted for. Also, the Bayesian updating used to study the strain rate scaling laws is briefly depicted. The following Section describes the macroscopic continuum model whose constitutive properties are obtained from nanoscale model. Numerical results will be presented before we conclude with a discussion in Section 4.

2 Nanoscale model

2.1 Model system and simulations

The material is described by a united atom (UA) model using the DREIDING force field [11] with harmonic covalent potential functions and the truncated Lennard-Jones (LJ) 6-12 for non-bonded van der Waals interactions whose parameters are adopted from [12]. The functional form and parameters are presented in Table 1.

The initial polymer structure was generated by using a Monte Carlo self-avoiding random walks algorithm as described by Binder [13]. A face-centered cubic (FCC) is used when generating initial configuration within a simulation box. Molecules were added to the lattice in a step-wise manner based on a method to make the appropriate selection of neighboring lattice sites. For each polymer chain, the first atom is added to an available site on the lattice. Then, the polymer chain is grown in certain directions on the bond angle and the density of the region where sites are not occupied in the probability context. LAMMPS [14] is employed to equilibrate the PE system through four sequential steps: (1) the PE structure was equilibrated for $10^{5}$ timesteps ($\Delta t=1$fs) at $500$K using a Nose-Hoover thermostat ($NVT$) [15, 16]; (2) a Nose-Hoover barostat ($NPT$) at the temperature of $500$K and the pressure of $1$atm condition was conducted for $5\times 10^{5}$ timesteps ($\Delta t=1$fs); (3) the structure was then cooled down to the desired temperature with a cooling rate of $0.4$K/ps followed by further $5\times 10^{5}$ timesteps ($\Delta t=1$fs) where the structure is in equilibrium. During the cooling process, the glass transition temperature ($T_{g}$) is determined as the intersection of two linear fitted lines to the volume versus temperature curve, see Figure 1(b). Three cooling rates $0.8$ K/ps, $0.4$ K/ps and $0.2$ K/ps are used herein to take the effect of cooling rate on the glass transition temperature ($T_{g}$) into consideration. As observed, volume-temperature plots and the resultant $T_{g}$ corresponding with various cooling rates are almost identical. It is shown that $T_{g}=300$K and density $\rho\approx 0.87\div 0.91$ g/cm3 are in good agreement with previous simulation and experiment results ($T_{g}=250K~{}\text{and}~{}\rho=0.95g/cm^{3}$ are experimentally measured value) on high density polyethylene (HDPE), see [17, 18, 19]. Furthermore, the influence of aging time on the stress-strain response was studied where the tensile stress-strain curves deformed at strain rate of $10^{10}~{}s^{-1}$ and temperature of $300$ K for three different polymer structures which are equilibrated by $500$ ps, $1000$ ps and $2000$ ps after the cooling process, respectively, are illustrated. As shown in Figure 1(c), in MD simulations when polymer systems are equilibrated long enough, the ageing time insignificantly influence on the stress-strain response as the curves are nearly the same for the initial stages. Deformation simulations will be described in the sequel.

2.2 Deformation simulations

In order to study the yield behavior of PE, the PE system was loaded in uni- and biaxially tensile/compressive strains at constant strain rates along the deformed directions. The pressure on the remaining two (uniaxial strain) or one (biaxial strain) lateral surfaces is maintained at $1$atm under NPT dynamics. As proposed previously [20, 21], the yield stress was taken as the maximum of stress-strain responses.

The Young’s modulus obtained from uniaxial tension at room temperature ($300$K) is $1.32$GPa, and the Poisson’s ratio is $0.32$, see Figure 2. These results are in good agreement with experimental results: Young’s modulus $E=1.38$GPa (obtained from testing method ASTM D368) and Poisson’s ratio $\nu=0.3$ [22, 23]. Note that the mechanical properties are averaged for three different initial PE structures to take entropic effects into account as suggested by [18]. Furthermore, the quasi-static tensile stress-strain response is simulated by using MD simulations as proposed by Capaldi et al. [17]. The system was uniaxially stretched at a constant strain rate of $10^{9}~{}s^{-1}$ for $1000$ steps followed by equilibration for $10000$ steps ($\Delta t=1$fs) with the axial dimension kept fixed to stabilize the energy in the system. This process is iterated until the desired strain is obtained. It is shown in Figure 2(c) that the quasi-static tensile yield stress ($\bullet$) is in a good agreement with experimental result reported by [23].

As the chain entanglement evolution is considered as important information that affects the deformation mechanisms of polymer, we have studied the chain entanglement evolution by using the geometric technique presented by Yashiro et al. [24]. As illustrated in Fig. 3(a), the interior angle $\theta$ between two vectors, i.e. one vector that is drawn from atom $i$ (A) to atom $(i-10)$ (B) and the other one that is drawn from atom $i$ (A) to atom $(i+10)$ (C), is measured. An example histogram of the distribution of the angles is shown in Fig. 3(b). The atoms, at which the angle $\theta$ is less than $90^{0}$, are classified as entangled or flexion nodes as indicated by [18]. Furthermore, the evolution of the entanglement parameter, which is obtained by dividing the number of atoms classified as entangled by the total number of applicable atoms, as a function of strain is plotted in Fig. 3(c). As can be seen, the entanglement parameter, which represents the percent of entangled atoms within the system, is nearly constant for the initial stages of deformation. At lager deformation ($\varepsilon\approx 0.5$), the entanglement parameter decreases nearly linearly with an increase in strain. These results are in good agreement with previous results reported by [18].

2.3 Evaluation of yield stress in multiaxial stress states

The principal stress components $\sigma_{i}$ are extracted from biaxially tensile and compressive deformations as proposed by [3]. Note that the stresses obtained from biaxial loadings have to be plotted versus the equivalent strain $\bm{\varepsilon}_{e}$ defined by [25]:

[TABLE]

where $\varepsilon_{ii}~{}\text{and}~{}\gamma_{ij},i,j=1,...,3$ are three normal and shear components of the strain tensor. The equivalent strain rate applied to the PE system is provided by Equation 1.

[TABLE]

Figure 4 shows the stresses in $x$ and $y$ directions versus the equivalent strain computed by Equation (1) and the corresponding yield peaks for biaxial tension and compression. The yield peaks ($\bullet$) occurring at the same equivalent strain are evaluated as maximum stresses.

2.4 Temperature dependence of elastic moduli

In order to predict the dependence of the Young’s modulus on the temperature the following Williams, Landel and Ferry (WLF) model [26] is used:

[TABLE]

where $T^{ref}$ is the reference temperature, $C_{1}$ and $C_{2}$ are adjustable WLF constants which are calibrated to fit the modulus data. Figure 5 shows that the Young’s moduli obtained from MD simulations are well explained by the WLF model.

2.5 Temperature dependence of yield stresses

Many studies have tried to account for the temperature and strain rate dependence of the yield behavior of polymers. The logarithm law [27], suggesting the stress-activated jumps of molecular segments results in yielding, is a good candidate to study the dependence of the yield stress on the temperature:

[TABLE]

where $\Delta H$ and $V^{*}$ are the respective activation energy and the activation volume; $\dot{\gamma}$ is the deformation rate, $\dot{\gamma}_{0}$ is a constant ($\dot{\gamma}_{0}\gg\dot{\gamma}$) [1]. The temperature dependent yielding law in Equation (4) can be approximately substituted by a linear fit (yield stress is considered as a linear function of the temperature) that is used hereafter. Cook et al. [2] also reported that the laws used to account for the dependence of yield behavior on the temperature for polymers are mostly linear.

Figure 6 shows the curves fitting the yield stress versus temperature for different strain rates are parallel. It means that the slope of the linear fits is nearly rate independent (even at quasi-static conditions). This has also been observed by Rottler et al. [28]. Hence, we are able to predict the yield behavior for different temperatures at quasi-static conditions, if the quasi-static yield stress value at any temperature is provided. Consequently, given the tensile yield stress ($25.92$MPa) at $T=300$K ($\bigstar$) obtained from MD simulations for quasi-static conditions, see Figure 2(c), the temperature dependent yielding law can be constructed at quasi-static strain rates (dash dot black line). A good agreement between this quasi-static linear fit and the experimental yield stresses extracted from [23] (with strain rate of $2~{}\text{min}^{-1}$) for different temperatures is quite clear. The predicted quasi-static tensile yield stress is related to the tensile yield stress obtained from MD simulations at the temperature of $300$K is expressed by

[TABLE]

The compressive and shear yield stresses as a function of the temperature for different strain rates are illustrated in Figure 7. The rate dependent the compressive and shear yielding laws obtained from MD simulations also show a parallel behavior. Furthermore, the tensile, compressive and shear laws (fitted lines) at different strain rates approximately change with the same rate suggesting the use of the same scaling factor $\gamma$ to predict the quasi-static compressive and shear yield stresses at the temperature of $300$K. The predicted quasi-static law for compression agrees well to experimental results obtained from ASTM tesing method [29, 22], as depicted in Figure 7.

[TABLE]

2.6 Strain rate dependence of the yield stress based on Bayesian approach

2.6.1 Bayesian updating

In this article, Bayesian approach is employed to calibrate the parameters of the plasticity constitutive models from the yield stress data obtained from MD simulations and existing experiments. The advantage of this method is that the naturally uncertain properties of the constitutive parameters existing in the multiscale model for polymers are taken into account. In the Bayes’ theorem the random variable $\bm{\theta}$ is expressed by a prior distribution $p(\bm{\theta})$. The uncertain parameters being estimated are then directly considered in the model evidence.

[TABLE]

with $\bm{\theta},~{}\bm{z}$ being the vector of model parameters and the vector of observations. For parameter identification purpose, the denominator $p(\bm{z})$ can be ignored and the posterior $p(\bm{\theta}|\bm{z})$is proportionally expressed by a combination of likelihood $p(\bm{z}|\bm{\theta})$ and the prior $p(\bm{\theta})$ as follows:

[TABLE]

Subsequently, the parameters yielding the maximum a posterior (MAP) probability of the parameters given the data is identified by:

[TABLE]

2.6.2 Scaling law constructed based on Bayesian approach

The strain rate in engineering practice is much lower than the one in MD simulations. Consequently, the respective yield stresses obtained from MD simulations and experiment can differ significantly. Hence, a scaling law is needed to upscale the yield behavior from nanoscale to macroscale. The good agreement between the predicted quasi-static yield stress ($\bullet$) and the experimental results at $T=300$K in Figure 2(c) implies that it is possible to rescale the high strain rate involved in MD simulations to macroscopic significant strain rate. Using the above-mentioned Bayesian approach, we can identify parameters of the strain rate dependent law when yield points at different strain rates obtained at molecular and continuum levels are determined. As suggested by earlier researchers [28, 30, 31], the dependence of the tensile yield stress $\sigma_{t}$ on the strain rate $\dot{\varepsilon}$ can be described by a logarithm or a power law form. In this article an exponential dependence of the yield stress $\sigma_{t}$ on the strain rate $\dot{\varepsilon}$ is adopted.

[TABLE]

where $\theta_{i},~{}i=1,..,4$ are constitutive parameters calibrated and updated on data points obtained from MD simulations and experimental data.

Since the temperature and strain rate are not correlated with respect to (w.r.t.) the yield stress as reported in [12], a linear transformation of the strain rate dependent yielding law (the numerical fit) for different temperatures is proposed. Richeton et al. [30] also suggested that the fitted law can be linearly transformed in vertical and horizontal directions when considering the effect of temperature and strain rate, respectively. Rate dependence of the yielding law is studied for different temperatures at the nanoscale model. Interestingly, as can be seen in Figure 8, fitted curves to the yield points at different temperatures are parallel. This supports the assumption that the rate dependent yielding law has also parallel behavior even at low strain rates. It means that under this assumption the rate dependence of yield stress for different temperatures can be predicted on a large ranges of strain rate (from low rate in practical application to high rate involved in MD simulations). As shown in Figure 9, predictions for the rate dependence of the yield stress at $T=250$K and $T=200$K are possible. The functional form in Equation (10) provides a good fit to the data. Good agreement of the predicted law with quasi-static yield stresses at $T=250$K and $T=200$K is observed. Hence, the model parameters can be identified in the case of limited experimental data from a Bayesian perspective. Note that the predicted quasi-static tensile yield stresses at $T=250$K ($\blacktriangledown$) and $T=200$K ($\blacklozenge$) correspond to the values illustrated by the same symbols ($\blacktriangledown$ and $\blacklozenge$) on the quasi-static fitted curve (black dash dot curve) in Figure 6, respectively.

3 Macroscopic Continum model

3.1 Definition of yield surface

The yield surface requires five constitutive properties: the uni- and biaxial tensile, compressive, and shear yield strength obtained from MD simulations. It is then constructed up to four Drucker-Prager-cones as suggested by Vogler et al. [4]:

[TABLE]

where $q=\sqrt{3J_{2}}=\sqrt{\frac{3}{2}\mathbf{s:s}}$ is the von-Mises equivalent stress; $p=-\frac{1}{3}I_{1}$ is the hydrostatic pressure, with ${I_{1}}=tr(\bm{\sigma})$ being the first stress invariant, and ${J_{2}}=\frac{1}{2}{\bm{s}}:{\bm{s}}$ being the second invariant of the deviatoric stress tensor $\bm{s}$ while $\varepsilon^{p}_{e}$ is the equivalent plastic strain. The parameter $\beta$ can be expressed in terms of the equivalent plastic strain [4] as

[TABLE]

where the parameters $\beta(\varepsilon^{p}_{e})$ is extracted from the hardening uniaxial tensile ($\sigma_{t}$), uniaxial compressive ($\sigma_{c}$), shear ($\sigma_{s}$), biaxial tensile ($\sigma_{bt}$), and biaxial compressive ($\sigma_{bc}$) curves which are obtained from MD simulations for corresponding stress states. The piecewise linear yield surface (PLYS) is illustrated in Figure 10.

In the proposed model, a non-associated flow rule is used to ensure the consistency with tensile test, leading to the plastic potential suggested by [10]:

[TABLE]

where $\alpha$ is the flow parameter accounting for the change in material volume at yielding:

[TABLE]

where $\nu_{p}$ denotes the plastic Poisson’s ratio obtained from the MD simulations under uniaxial tension and

[TABLE]

The increment of plastic deformation is given by:

[TABLE]

where $\Delta\lambda$ is the plastic multiplier, commonly updated via the return mapping algorithm under the Kuhn-Tucker consistency conditions. A more efficient approach based on Chen Mangasarian replacement functions which avoids a return mapping has been proposed by [32, 33]; $\frac{\partial g}{\partial\bm{\sigma}}$ represents the direction of plastic flow with $g$ being the plastic potential given in Equation (13). The equivalent plastic strain is given by [34]:

[TABLE]

with $k=\frac{1}{1+2\nu_{p}^{2}}$.

3.2 Thermo-plastic hardening

The constitutive model is defined by uniaxial tension and compression, biaxial tension and compression and shear yield strengths. Thus, the hardening will be formulated to update these yield strengths. Commonly to other plasticity models, the hardening formulation depend on the equivalent plastic strain as follows:

[TABLE]

We can directly extract stress and strain values from the uni- and biaxial tension and compression and shear from MD simulations. Then, the hardening laws were inserted into the material model in terms of table of values. Note that input data are presented in terms of plastic strain by decomposing the total strain increment by the elastic component as: $\Delta\varepsilon^{p}=\Delta\varepsilon-\Delta\varepsilon^{el}$ [10]. In each iteration the table lookups will provide the plastic strains $(\varepsilon^{p}_{e})$ and corresponding yield stresses $(\sigma_{y})$ as inputs. Subsequently, the tangents $\left(\frac{\partial\sigma_{y}}{\partial\varepsilon^{p}_{e}}\right)$ with respect to the plastic strain will be computed. These stress-plastic strain curves are then scaled to determine the hardening laws. In order to study the temperature dependent yield strength, the linear law fitted on data obtained from MD simulations is employed to scale the yield stress w.r.t. the hardening curve at the reference temperature as follows:

[TABLE]

with $\sigma_{n}$ and $\sigma^{ref}_{n}$ being the predicted yield stresses at the desired $T$ and reference $T^{ref}$ temperatures, respectively. The material constant $\beta_{n}$ is selected so that the yield stresses at the temperature $T$ are scaled back to the stresses’ value at the reference temperature $T^{ref}$. An overview illustrating the algorithm that is applied to implement the PLYS constitutive model is shown in Table 2.

3.3 Yield surface at different temperatures

We perform multiaxial deformations (uni- and biaxial and shear loads) to obtain yield points at two different temperatures and the equivalent strain rate in Equation (1) is set as $\dot{\varepsilon}_{e}=1\times 10^{10}~{}s^{-1}$. The PLYS characterized by Equation (11) was adopted to fit yield points data in four Drucker-Prager cones as mentioned in Equation (12). As can be seen in Figure 11, the yield points are well described by the PLYS criterion.

3.4 Yield surface at different strain rates

Based on the scaling law proposed in Section 2.6, given any known (predicted) yield stress of a specific load case, the entire multisurface yield functions can be isotropically scaled to quasi-static rates by assuming the scaling value is similar for general deformations. For example, a prediction for the entire quasi-static multisurface yield functions at $T=300$K is obtained in Figure 12(a). The uniaxially quasi-static tensile ($25.92$MPa) and compressive ($37.6$MPa) yield stresses are validated with experimental results. As observed in Figure 12(a) and Table 3, good agreement between numerical results and experimental results is observed. Furthermore, the entire yield surface at any desired strain rate can also be predicted based on the law shown in Figure 9 using Equations (5 + 6). Also, the entire yield surfaces at $T=250$K for different strain rates are obtained from MD simulations and the one at quasi-static rates is predicted using the same scaling law as illustrated in Figure 12(b).

The above-described elasto-plastic model is used to predict the thermoplastic behavior at (1) nanoscale and (2) the strain rate, which is rescaled from molecular to continuum levels through the constitutive law.

The presented constitutive model in the aforementioned section was implemented as material parameters into ABAQUS to predict the macroscopic stress-strain responses. Comparison between the responses obtained from MD simulations and from the continuum model for different stress states in Figure 13 shows a good agreement. The temperature dependence of the uniaxial and biaxial tensile, compressive and shear stress-strain responses is also illustrated in Figure 14 and good agreement between the responses obtained from the continuum model and those from MD simulations for different temperatures is observed. Furthermore, the consistency of the stress-strain responses for different stress states (e.g. the unequally biaxial tension and compression with $\dot{\varepsilon}_{y}=2\dot{\varepsilon}_{x}$ yield stresses) between MD simulations and the continuum model could be expected.

Figure 15 compares the predicted yield strength ($\sigma^{quasi}_{t}=27.65$ MPa) with the experimental one reported in [23] ($\sigma^{exptl}_{t}=29.3$ MPa) at strain rate of $2~{}\text{min}^{-1}$. Also, the continuum model accurately predicts the compressive yield stress at room temperature, see Table 3. This proves that the continuum model can be used to accurately predict the yielding occurring at low strain rate.

4 Conclusions

A hierarchical multiscale model was developed to study the thermo/visco-plastic behavior of the PE. At first, the PLYS and the temperature and strain rate dependent yielding laws were constructed where the constitutive parameters are calibrated from data (yield points for multiaxial stress states) obtained from MD simulations. Then, the scaling law for the entire yield surface was proposed based on the quasi-static tensile simulations at nanoscale. The yield behavior was upscaled to macroscopic level through an efficient continuum model. The consistency of the results demonstrates that the macroscopic continuum model accurately predicts the behavior achieved from MD simulations.

In addition, validation shows that the tensile and compressive yield stresses are accurately predicted at quasi-static rates by the proposed multiscale multisurface model despite the loss of ad-hoc experimentation. Hence, we believe that this study will open a new door for the design of polymer materials through multiscale simulations, leading to priori predictions of yield behavior of polymers.

5 Acknowledgements

We gratefully acknowledge the support by ERC COMBAT project (project number 615132).