Probabilistic estimation of the constitutive parameters of polymers

The Mulliken-Boyce constitutive model predicts the dynamic response of crystalline polymers as a function of strain rate and temperature. This paper describes the Mulliken-Boyce model-based estimation of the constitutive parameters in a Bayesian probabilistic framework. Experimental data from dynamic mechanical analysis and dynamic compression of PVC samples over a wide range of strain rates are analyzed. Both experimental uncertainty and natural variations in the material properties are simultaneously considered as independent and joint distributions; the posterior probability distributions are shown and compared with prior estimates of the material constitutive parameters. Additionally, particular statistical distributions are shown to be effective at capturing the rate and temperature dependence of internal phase transitions in DMA data.


Introduction
Constitutive relationships are typically developed to capture the strain-, temperature-, and rate-dependent response of various materials. This paper examines the rate and temperature dependence of dynamic mechanical analysis (DMA) and dynamic compression data for polyvinyl chloride (PVC). The two principle objectives are (1) propose a new parameterization of the phase transitions in DMA data, and (2) establish sensitivity to uncertainty in the estimation of constitutive model parameters of the Mulliken-Boyce constitutive model for crystalline polymers [1] using Bayesian inference.

Mulliken-Boyce model
The Mulliken-Boyce (M-B) model is a two phase (α and β) viscoelastic-viscoplastic model, shown schematically as Maxwell-Weichert elements in Figure 1, that incorporates both the polymer network stress (B) and the chain stresses (A) due to the polymerization. The M-B model has been shown to capture rate-dependent behavior of polymers [2,3] very accurately, especially the rate-dependent yield strength and post-yield response.

Governing equations
The Mulliken-Boyce model proposes two contributions to the total stress in the material. The intermolecular stress (T A ) resists chain segment rotation and has two principle rate-activated phases, α and β. Both phases are modeled as viscoelastic-viscoplastic processes ( Figure 1). The network back stress (T B ) along the polymer chain resists alignment and is modeled as a Langevin stress [4]. The total stress in the system is therefore Using a 1-D approximation, the total stress is σ total = σ A,α + σ A,β + σ B (2) and the strain is It follows that the strain rate iṡ The 1-D expression for the Langevin (network) back stress is where λ = U 2 + 2U −1 /3 1 /2 is the chain stretch parameter in 1-D, U = exp(ε), and L −1 inverse of the Langevin function, A Padé approximation [5,6] is used to calculate the inverse Langevin function, Under the assumption of uniaxial response, the linear system of governing ordinary differential equations becomeṡ This is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.  It is important to note that the strain is carried forward as a time-dependent parameter in contrast with the standard practice of assuming constant strain rate. The remaining parameters are explained in Table 1. This addresses the practical issues of the early time response of polymers to incident stress waves where the system has not yet reached dynamic equilibrium and non-ideal experiments where the rate evolves continuously (i.e., never reaches equilibrium). Additionally, this approach enables the use of implicitly rate-dependent parameters (currently being explored as an improvement to the solution).

Solution of ODE
The nonlinear system of equations in Eq. (8) is solved using Matlab ODE solvers. The multistep ode113 function, based on the Adams-Bashforth-Moulton PECE solver [7], is used due to the computationally expensive nonlinear relationships. An example predicted response for PVC is shown in Figure 2 at a relatively low strain rate (ε = 10 s −1 ) for three temperatures.

Bayesian inference of parameters
Constitutive parameter estimation is an inverse problem wherein experimental data (observations) are used to infer the underlying parameter value(s) [8]. Bayesian statistics provide an intuitive, fundamental construct for the consideration of uncertainty in modeling, independent of its source (e.g., random experimental error vs. process-driven uncertainties). In this Bayesian framework, the conditional probability of the unknown parameters (x) based on the observed data (d) is the posterior probability distribution function (p) and is given for normally distributed data by where m(x) is the model output based on the unknowns, x 0 is the prior value for the unknown parameters, Γ 0 is the covariance matrix, Γ d is the covariance matrix of the data, and N is a normalization constant. Maximizing Eq. (9) is accomplished using a maximum a posteriori (MAP) estimator.
In the absence of available information, Bayes' postulate states that a uniform prior probability should be assumed [9]. In the case of continuous variables, a uniform pdf over all x (called a diffuse, non-informative, or Dirichlet prior condition) is used, i.e., where |x| < (2ε x ) −1 and is ε −1 x sufficiently large to span the plausible range of x.
Two different experimental data sets of PVC are analyzed, DMA and dynamic compression. The underlying parameters estimated are experiment-specific; these are reviewed in the next two sections.

Sample material and geometry
PVC was chosen for study due to the availability of relevant data for comparison [10,11]. Impact resistant PVC (Type II) was machined from extruded 25.4 mm diameter rods into samples that measured 60 mm long, 12.5 mm wide, and 3.2 mm thick. These samples were tested in a dual cantilever configuration in a TA Instruments Q800 [12,13] at frequencies of 1, 10, and 100 Hz and a temperature range of −100 • C to 190 • C. The displacement was held constant at 15 µm for this analysis.

Experimental DMA data and processing
Typical DMA data for the PVC is shown in Figure 3 at the measured frequencies of 1, 10, and 100 Hz. The frequency is converted to strain rate (ε) equivalence relationship . [1], where ω is the angular frequency (in rad/s), d 0 is the amplitude of the displacement, and l g is the specimen gage length. The corresponding time-temperature superposition for strain is typically included in the analysis using the WLF equation [14], where T 0 and ω 0 are reference quantities and C 1 and C 2 are constants. An alternative expression for he DMA shift [15,16] is also considered, which parameterizes the shift using only two parameters, T 0 and the constant A.

Phase transition statistical model
In the Mulliken-Boyce model, each of the two phases (α and β) are due to unique physical mechanisms. In Figure 3, the loss modulus show two distinct transitions that can be attributed to the phase: the α transition (i.e., the glass transition) occurs at T ≈ 350 K and the β transition occurs at T ≈ 200 K to 250 K. Since both the α and β transitions are spread over a finite temperature range, it is appropriate to consider each parameter as a statistical variable with underlying pdf's. The form of these distributions for α and β processes are modeled differently. This is supported due to the statistical mechanics of the underlying transition [17], as discussed in the following sections. Additionally, the two phases' properties, such as transition temperatures and activation energy, are assumed to be independent (not causally related).

Alpha phase
The loss modulus of the glass transition (α phase) gives an indication of the irreversible conversion of mechanical energy into heat by the breaking of inter-chain bonds [17]. In a statistical interpretation, the loss modulus probability distribution p E α is the likelihood of the polymer bond being broken at a given temperature. The loss modulus of the α phase (E α ) is further assumed to be proportional to the temperature gradient of the storage modulus α phase The associated cumulative distribution function, χ E α of the bond breaking distribution is given by which can be interpreted as the accumulated loss of polymer strength (i.e., storage modulus) due to the cumulative effect of bonds as the system heats through the glass transition. This is taken to be proportional to the storage modulus probability distribution, p E α

EPJ Web of Conferences
Mahieux [18] proposes a Weibull distribution [19] for the modulus of a polymer with respect to temperature, W(T ) due to the cascading interaction of ruptured bonds as it undergoes a glass transition. The probability distribution of the loss modulus, E α then assumes the form where ∆E α is a modulus contribution by the α phase, T α is the transition temperature, and is the distribution shape parameter. It follows that the storage modulus distribution is given by where A α is a constant. The reduced storage modulus with increasing temperature can then be regarded as a progressive, fatigue-like accumulation of damage in the polymer network that reduces the polymer's strength (due to the reduced number of available bonds).

Beta phase
For the β phase, the same reasoning applies to the pdfto-cdf relationship between the loss and storage modulus. However, the form of the underlying distribution is different since the β phase is a normal distribution [17], i.e., where ∆E β is the modulus contribution of the β phase, is T β the transition temperature, and α T β is the distribution width. The storage modulus of the β phase is then

Storage and loss moduli
The overall model for the loss modulus is the superposition of the contributions from each of the two phases with an additional normally distributed residual (or background) loss modulus, i.e., The storage modulus is similarly constructed,

Model-based estimation
The unknown parameters to model the DMA moduli distributions are given in Table 2. The corresponding vector of unknowns is The best unknown parameter distributions are initially found using a stochastic multistart global optimization routine in Matlab [20]. These points then become the initial guesses for a deterministic MAP estimator to refine the estimates further. The results of fitting are shown in Figure 4. It is important to note that while the best observed fit (with maximum posterior probability) is not necessarily a global minimum, the large number of estimates minimizes the likelihood that the global minimum was missed.
A Monte Carlo approach [19] is implemented to characterize the posterior distributions of the estimates themselves and to also characterize the sensitivity of the model to variations in the prior distributions. width of the storage modulus, with σ T β = 5 K. Figure 6 shows the posterior distributions for the first six elements of x. The diagonal elements show the posterior distributions for each of the parameters while the off-diagonal terms are indicative of the statistical covariance between the estimation variables. For example, the posterior distributions of ∆E α (x 1 ) and (x 4 ) are normally distributed while other elements of x appear bimodal. The covariance also indicates that ∆E α and ∆E β are linearly related, which would be expected since both parameters are assumed to contribute to the storage and loss moduli via superposition. The underlying mechanisms for the dependencies between all of the estimation parameters are currently being explored.

Discussion of DMA results
The ability to generate and estimate the statistical distributions for the phase contributions to the modulus based on polymer chain dynamics is encouraging. It should be emphasized, however, that the estimates are based on a single sample, i.e., a single realization of the material. Further work remains to determine the distribution uncertainties as a function of intrinsic material variability or other external factors, such as processing.
The analysis in this paper assumes uniaxial behavior. If symmetry is assumed, then only a single β mode would be expected. However, the DMA method is based on a finite cantilevered beam in bending; this creates the possibility that there are in fact two slightly different natural axes due to the different levels of stress in the principle (orthogonal) axes of the beam. In this case, it is relevant to entertain the idea of multiple β contributions with closely spaced but independent parameters. The skeletal modes would then be considered as β modes, and by superposition the loss and storage moduli would be and (25) respectively. The extra four parameters enable a more accurate fit of the data and will be examined in the future.

Dynamic compression of PVC
The other experimental method considered is dynamic compression experiments of PVC. The methods and data are reviewed in Jordan et al. [21]. The unknown model parameters in the Mulliken-Boyce model estimated from the dynamic compression experiment data are The DMA storage modulus data is used to generate the elastic modulus for the model. Other parameters are based on results in the literature [1,10]. Iterative estimation from multiple data sets and different strain rates were performed using a global optimization scheme with different initial guesses. Figure 7 shows results from the estimation of PVC compression data at one strain rate (ε = 0.01s −1 ) Consistent with previous efforts [3], low rate experiments yielded better estimates of the α phase parameters whereas the β phase was only active in high rate experiments and cannot be accurately estimated with low rate data. The posterior distribution of the α phase parameters is shown in Figure 8. As with the DMA estimates, the posterior distributions vary significantly and the underlying causes are being explored. Additionally, the strong covariances observed in the off-diagonal terms in Figure 8 imply that additional consideration to covariance between parameters may be warranted.

Conclusions
A Bayesian framework was implemented for the Mulliken  developed. A Weibull distribution was used to fit the α phase parameters whereas a normal distribution was used to model the β phase. The constitutive parameters of the Mulliken-Boyce model in dynamic compression experiments were estimated. Finally, the posterior distributions and model sensitivity to the underlying distributions of both experiments were characterized using a Monte Carlo approach.