Predicting embrittlement of polymer glasses using a hydrostatic stress criterion

In this study, the aging-induced embrittlement of three polymer glasses is investigated using a previously developed hybrid experimental–numerical method. The evolution of yield stress of unnotched tensile bars upon aging is coupled to the evolution of embrittlement of notched tensile bars using a numerical model combined with a critical hydrostatic stress criterion that determines the onset of failure. The time-to-embrittlement of notched tensile bars with a different notch geometry is predicted and in good agreement with the experimentally determined value. Next to that, the approach is extended to three polysulfone polymers, and it is shown that the value of the critical hydrostatic stress correlates well with the polymers entanglement density: : polymers with a denser entangled network display higher values, that is, a higher resistance against incipient cavitation. © 2019 The Authors. Journal of Applied Polymer Science published by Wiley Periodicals, Inc. J. Appl. Polym. Sci. 2019, 136, 47373.


INTRODUCTION
Amorphous polysulfones (PSUs) are rigid and tough materials that are therefore widely employed in engineering applications for their outstanding properties. However, in time, the toughness of polymers can be greatly reduced as a result of physical aging, and a transition from a ductile to a brittle failure mode, that is, embrittlement, may be observed. 1 The understanding of this transition is of great importance for polymer components used in load-bearing applications where ductile failure is usually preferred. Furthermore, despite of similarities between the molecular architecture of PSUs, there are marked differences in their impact response. [2][3][4] A method eminently suitable to quantify differences in the impact response was demonstrated by Engels et al. 5 In this method, the evolution of yield stress upon aging of unnotched tensile bars is linked directly to the time-dependent embrittlement of notched tensile bars. The presence of a notch results in severe localization of deformation with increasing strain and the resulting gradients in volumetric strains lead to a stronger build-up of a hydrostatic stress underneath the notch. At some point, the hydrostatic stress reaches a critical level, and a void is created to release the stress. These voids initiate a craze which ultimately triggers brittle failure. [5][6][7][8][9] The build-up of hydrostatic stress is a local phenomenon and is therefore not directly accessible in experiments. Narisawa and coworkers were the first to identify this local hydrostatic stress for several polymers using a numerical analysis employing the slip-line theory. [10][11][12] Studies of van Melick et al. 9,13 convincingly demonstrated that the build-up of hydrostatic stress can also be successfully analyzed in finite element (FE) simulations using a constitutive model that captures a polymers intrinsic deformation response. Similar results were reported by Gearing and Anand. 14 It is well-known that polymer glasses are typically not in thermodynamic equilibrium and as a result will display a drive toward equilibrium, that is, physical aging. This leads to a gradual change in mechanical properties over time, as for example, an increase in elastic modulus or yield stress. 15 The typical increase in yield stress upon physical aging induces an increase in strain softening, leading to a stronger tendency to plastic strain localization that ultimately leads to embrittlement. 16,17 Engels et al. 5 showed that by monitoring both the yield stress evolution and embrittlement upon annealing that the transition from a ductile to a brittle failure mode corresponds to a critical thermodynamic state of the material, reflected in a critical value of the evolving yield stress. Furthermore, they showed, using FE simulations, that the increase in yield stress upon annealing directly translates to an increase in the maximum value of the hydrostatic stress behind the notch during impact. This implies that the experimentally obtained critical yield stress can be directly related to a critical hydrostatic stress. Interestingly, the critical hydrostatic stress was shown to be material specific and was suggested to be related to the entanglement density 9, 17,18 ; where a high entanglement density is observed to give a higher resistance against voiding, that is, a higher value of the critical hydrostatic stress.
Several constitutive models have been reported in literature to numerically simulate localization phenomena. Most are based on the work of Haward and Thackray, who were the first to recognize two separate contributions to the stress response of a polymeric material: a viscous contribution of the intermolecular interactions and an entropic contribution related to the orienting entangled molecular network. 19 Examples of such constitutive models are the BPA model, 20,21 the Oxford glass rubber model, [22][23][24] and the Eindhoven glassy polymer (EGP) model. [25][26][27][28] Despite of differences in their approaches, all these models demonstrate the significance of the intrinsic material response on plastic strain localization. [29][30][31][32] In the most recent version, the EGP model is able to account for rate-and temperature-dependent strain hardening, achieved by a means of a viscous contribution to the strain hardening stress besides an entropic contribution. 28 This viscous contribution expresses itself as a deformation dependence of the flow stress, as proposed earlier by Wendlandt et al. 33,34 The influence of physical aging is captured in the EGP model by a state parameter S A , which uniquely describes the thermodynamic state of the material and includes strain-induced mechanical rejuvenation. 26 The goal of the present study is to use the framework presented by Engels et al. 5 to evaluate aging-induced embrittlement and the predictability thereof in combination with the latest version of the EGP model presented by Senden et al. 28 The PSUs of interest are three commercially available systems: polyphenylsulfone (PPSU), polyethersulfone (PESU), and PSU. EGP model parameters for each material were derived from the intrinsic material response, measured in uniaxial compression tests at various strain rates and temperatures. The annealing-induced embrittlement was investigated in uniaxial tensile experiments on notched samples. The advantage of uniaxial tensile testing over Charpy or Izod impact testing is that the boundary conditions are better defined, which enable more straightforward FE modeling. FE simulations of the notched tensile test were used to simulate the build-up of the local hydrostatic stress, using the experimentally obtained critical yield stress, reflected in the value of the state parameter, as a measure for the thermodynamic state of the material. The obtained critical hydrostatic stresses are compared with literature data on other amorphous polymers and are in full agreement with the previously reported correlation with entanglement density. Furthermore, using the hydrostatic stress as a criterion, the time-to-embrittlement of notched tensile bars with a different notch geometry could be predicted.

Materials
Materials used are Radel PPSU, Udel PSU, and Veradel PESU kindly supplied as granules by Solvay Specialty Polymers (Alpharetta, GA). Table I gives an overview of the different grades and corresponding molecular weights. Extruded rod of PPSU and PSU was obtained via Solvay or commercial partners.
Tensile bars were injection molded according to ASTM D 638 Type I specifications (dimensions: 100 mm gauge length, 3 mm thickness, and 13 mm width) using an Arburg Allrounder 320S 500-150, using a mold temperatures of 150 C. Initially, this mold temperature was used for PESU but subsequent impact tests indicated that the samples displayed a brittle failure mode. Consequently, the mold temperature was reduced to 55 C to obtain samples with a lower yield stress and a ductile failure mode. Processing conditions were kept the same for different grades of the same material. Notched tensile bars were obtained by a milling operation at room temperature. Three different notch radii were used: 0.25 and 0.4 for PPSU; 0.8 mm for PSU and PESU. All notches have an included angle of 35 and a depth of 3.5 mm.
For uniaxial compression samples, plates of 1 cm thick were compression molded at temperatures of 350 C (PSU) to 370 C (PPSU and PESU). Granules were allowed to melt in a mold covered with aluminum foil at atmospheric pressure. Subsequently, the material was compressed up to 100 kN for 2 min and the mold was cooled down in a cold press at 20 C, applying little pressure. As the molding procedure might result in a surface layer with different mechanical properties, these sections are removed by reducing the thickness of the plate symmetrically to 7 mm by a milling operation. Cylindrical samples measuring 6 × 6 mm (diameter × height) were subsequently machined from the plates on a turning lathe. To determine the pressure dependency parameter compression samples of 4 × 4 mm, as well as small tensile bars, were milled from the center of another plate to assure both sample geometries experienced the exact same thermal history.
Annealing treatments were performed using hot air circulating ovens at various annealing temperatures and annealing times. Prior to testing, annealed samples were allowed to cool down to room temperature.

Experimental Methods
Uniaxial compression tests were performed on a Zwick 1475 tensile tester equipped with a 100 kN load cell at several true strain rates and temperatures (using an MTS temperature chamber). Friction between sample and steel plates was reduced by applying polytetrafluoroethylene (PTFE) tape (3M 4580) to the sample and PTFE spray (Griffon TF089) to the steel plates. True strain rate control was used, under the assumption of incompressibility, that is, a constant rate of true strain, and true stress and true strain  Figure 1.

Constitutive Model
The EGP model 26 is used in its multimode, 27 multiprocess 35 form with an extension to incorporate a viscous contribution to strain hardening. 28 For completeness, a summary of the governing equations is given here. In the model, the total Cauchy stress is additively split in a hydrostatic component and a deviatoric component: The hydrostatic component described with a constant bulk modulus κ and the volume change ratio J: The deviatoric component consists of an elastoviscoplastic driving stress and an elastic hardening stress: where the strain hardening stress is described by the Edwards-Vilgis model, for details see Refs. 31, 36 and 37. For the driving stress, thermorheological complex behavior is taken into account: two relaxation processes are assumed to work in parallel, an α-process (related to the primary glass transition, that is, mainchain segmental motion) and a β-process (related to the secondary glassy transition, that is, partial main-chain or side-chain motion), each with n and m individual modes respectively. A primary and secondary transition is typically observed in PSUs. 38 Furthermore, it is assumed that the stress contribution of each relaxation process can be described by an arbitrary number of discrete viscoelastic Maxwell modes, that is, using a spectrum of relaxation times. The expression for the driving stress then reads: withB d ex,i the deviatoric part of the elastic, isochoric left Cauchy-Green strain tensor and G x, j the shear modulus. Note that, subscript x refers to a specific process and i to a specific mode. Because of kinematic considerations the plastic rate of deformation tensor D p has to be specified. This is achieved with a non-Newtonian flow rule relating D p to the driving stress with a viscosity η x, i : The viscosity is function of the equivalent shear stress τ x and absolute temperature T and is described by an Eyring type of relation. Furthermore, it is extended to take hydrostatic pressure p into account: Here, η 0x,i is the initial viscosity, τ 0x is the characteristic shear stress, ΔH x is the activation enthalpy, R is the universal gas constant, and μ x is the pressure dependence parameter. The viscous contribution to the strain hardening response expresses itself by a deformation-dependent activation enthalpy and initial viscosity, using invariant functions of the total strain I rB À Á , withB the total isochoric elastic left Cauchy-Green strain tensor. 28 The characteristic shear stress τ 0x , the pressure p, and the equivalent shear stress τ x are given by: in which k B is the Boltzmann constant and V * x is the shear equivalent activation volume. Intrinsic strain softening is captured with the state parameter S x : The initial value S aα uniquely defines the thermodynamic state of the material. Note that the state parameter of the β-process is taken as zero. The softening function, R γ , which is a function of equivalent plastic strain γ p , captures the softening kinetics and is given by 26 : where r 0,x , r 1,x , and r 2,x are constants. The equivalent plastic strain γ p is coupled to the mode with the highest relaxation time, that is, the α-mode, where i = 1, and can be derived from the corresponding equivalent plastic strain rate according to: Material Characterization To obtain all model parameters for the EGP model several steps have to be performed, which will be shortly discussed below. Note that Points 1-4 are covered in a different publication, 39 and their determination will not be discussed here in detail. These are the parameters that mainly describe the large strain response, namely, the elastic strain hardening stress and the viscous contribution to strain hardening, which are both not affected by the thermodynamic state of the material. The steps to be performed are: 1. The pressure dependency parameter μ x , which captures the influence of hydrostatic pressure p, is determined from a series of compression and tensile tests on samples with similar thermal histories at several strain rates. The difference between the tensile and compressive yield stresses is in this case fully related to the difference in hydrostatic pressure. Using a pressure-modified Eyring expression, 40,41 the pressure dependency parameter can be determined. 2. The elastic hardening parameters, G r and α r , a parameter related to the limited extensibility of the network, as well as the bulk modulus κ, total shear modulus G, and total initial viscosity η 0 are determined from either a tensile test on a mechanically rejuvenated axisymmetric tensile bar or a compression test. 3. The distribution of elastic and viscous hardening is determined following Senden et al. 42 by performing compression tests on predeformed samples. 4. Subsequently, the parameters describing the viscous contribution to strain hardening can be determined following Senden et al. 28

The softening parameters are determined by fitting a modified
Carreau-Yasuda function to the postyield softening response using a single uniaxial compression test following Klompen et al. 26 6. The multimode relaxation spectrum is derived from the preyield response of a single compression test following van Breemen et al. 27 A spectrum of relaxation modes yields an improved description of the preyield response, while the postyield response remains unaffected. 27,36 To the obtained α-spectrum, a single β-mode is added.
The parameters in Points 5 and 6 are important for the preyield response and postyield softening, and thus strain localization which evidently is of importance for the type of simulations considered in this article. Furthermore, they are affected by the thermodynamic state of the material. First, the softening parameters (Point 4) are determined. To do so, the driving stress σ s is isolated by subtracting the hardening stress from the total stress. The resulting driving stress as shown in Figure 2(a) for PESU, can be split in a rejuvenated stress σ rej and an aging-dependent yield drop Δσ y . Assuming that at the yield point, the equivalent plastic strain γ p starts to evolve and that the softening function equals 1, we can write: The resulting softening characteristic is shown in Figure 2(b) for all three materials. Solids lines are fits with eq. (10), markers represent experimental data and a good fit is obtained. It should be noted that the softening characteristic is assumed not to depend on strain rate or temperature, which appears to be valid for the range of strain rates and temperatures used in this study. Furthermore, the parameters are assumed to apply to the β-process as well, since no data were available in this regime to determine the β-process softening characteristic.
Next, the multimode relaxation spectrum is determined from the preyield response of a single compression test (Point 5). It has been shown that a multimode spectrum improves the preyield description which is of utmost importance for uniaxial impact simulations. 27 The method used relies on the time-stress superposition principle, and more detailed information can be found in van Breemen et al. 27 The multimode EGP model consists of an arbitrary number of generalized, nonlinear Maxwell elements, that are assumed to work in parallel and are fully related to the α-process. The constitutive behavior of such a Maxwell element, in 1D, can be expressed by a Boltzmann integral in its relaxation form: where σ s (t) is the driving stress at time t, E is the relaxation modulus, and _ ε is the strain rate. The stress reduced time ψ and the pending stress reduced time ψ 0 can be calculated by integration of the stress shift factor a σ (σ), 36,43 which implies that the relaxation time of the Maxwell elements becomes shorter when a stress is applied. For a multimode spectrum, the relaxation modulus is expressed as: with λ i the relaxation time of mode i. Substitution of this equation into the previous yields: The multimode spectrum is thus characterized by a discrete spectrum of relaxation times. To obtain the relaxation spectrum, the driving stress up to yield is isolated from the total stress, measured in a single compression test, by subtracting the elastic hardening stress and the above integral is evaluated at every experimental time point for each separate relaxation time. The experimental curve is obtained at a temperature of 22 C, at a strain rate of 10 −3 s −1 , where no influence of the β-process is observed. The obtained spectrum is thus fully related to the α-process. The number of relaxation times necessary is somewhat arbitrary, but in general one relaxation time per decade of time is necessary to obtain an accurate description. A too low number of modes results in a nonsmooth relaxation curve while a too high number will only result in excessive computation times. The number of relaxation times, and thus the number of discrete modes, is 23, 17, and 19 for PPSU, PSU, and PESU, respectively. The obtained elastic relaxation moduli are converted to shear moduli G α,i , see Ref. 27, and the initial viscosities are subsequently calculated using η * 0α,i ¼ λ i × G α,i . The obtained spectrum is still related to an aged material and therefore the calculated viscosities have to be corrected for the current thermodynamic state of the material. This is done by equally shifting all viscosities along the time axis by a horizontal shift only (time-aging time superposition principle), implying that all the relaxation times are equally affected by physical aging 15,27,44 : To the obtained spectra, the single shear modulus and initial viscosity of the β-process, taken from Ref. 39 are added yielding a total number of modes of 24, 18, and 20 for PPSU, PSU, and PESU, respectively. The obtained spectra for all three polymers are listed in the Appendix.
In Figure 3, the influence of thermal history, achieved by severely annealing some of the compression samples, on both the relaxation spectra (left figures) and the corresponding intrinsic stress-strain response (right figures) of all polymers is shown. The rejuvenated responses, corresponding to the reference spectrum, are indicated with the down triangles in all figures and are obtained using a state parameter equal to zero, that is, no strain softening. To describe the stress-strain response of the material with a different thermomechanical state, the reference spectrum is shifted with the appropriate value of S a (solid lines in Figure 3, left). To show that this shift of the reference spectrum works well, the spectrum of the materials with a different thermomechanical state are also determined using the method describe above (dashed lines), and as can be seen these coincide indicating that time-stress superposition is an adequate approximation for the present purpose. With the shifted reference spectra, the corresponding stress-strain curves are calculated and compared with experimental data; results are shown on the right-hand side of Figure 3. With an increase in annealing temperature (annealing time is 24 h), both the yield stress and strain softening (the decrease in stress after yield) increase whereas the large strain response remains unaffected by aging. As can be seen, the obtained softening parameters describe the strain softening response accurately. With an increase in yield stress upon aging, the value of the state parameter increases likewise and the shift of the reference spectrum accurately describes the preyield response for PPSU and PESU. For PSU, the magnitude of the yield stress is captured accurately; however, for the severely aged sample, physical aging does not only result in an increase of the yield stress (and strain softening), but also the modulus is affected. With annealing, the modulus displays a stronger increase than observed in the other two materials. This mismatch becomes also apparent in the relaxation spectrum, where the initial modulus of the aged material (dashed line) is higher than the one of the shifted reference spectrum (solid line). To capture the increase in modulus in PSU, a so-called vertical shift, which is not uncommon, 15,44 would also be necessary to fit the experimental preyield response of the aged samples. It appeared that this can effectively be achieved by increasing the value of the first α-mode shear modulus (G α,1 ) from 325 to 520 MPa, to fit the experimental preyield data. For all other EGP simulations on PSU, this increase in modulus is checked and if needed accounted for by changing the value of (G α,1 ). For completeness, all model parameters for the three polymers can be found in the Appendix.
The obtained material parameters are further validated by performing uniaxial compression simulations at several strain rates and compare the results to experimental data. In Figure 4(a), the compressive stress-strain curves are shown, measured at room temperature. The EGP model accurately describes the uniaxial compression response of the three polymers. In comparison, PPSU displays the most strain hardening, also indicated by the value of the elastic hardening modulus G r which equals 8 MPa, followed by PSU (G r = 5.3 MPa) and PESU (G r = 4 MPa). Interestingly, the drop in stress after yield follows the opposite trend (see also Figure 2). The balance between strain hardening and strain softening is believed to be important for the toughness of the polymer: moderate strain softening and pronounced strain hardening are typical of a tough material, for example, polycarbonate, whereas pronounced strain softening and only moderate strain hardening are typical for a brittle material, for example, polystyrene. 17 Figure 4(b) shows the rate dependence of the yield stress measured in uniaxial compression at room temperature. Solid lines are EGP simulations. With an increase in strain rate, the typical increase of yield stress is observed and captured accurately by the EGP model.

Quantifying Aging-Induced Embrittlement
It is well-known that the crazing behavior of a polymer is controlled by the entanglement density. 45 By altering the network density of PS-PPE through variation of the PS/PPE ratio in the blend, van Melick et al. 9 showed that the critical hydrostatic stress increases with increasing network density. The hydrostatic stresses reported were obtained with a previous version of the EGP model. To further investigate this relation, the critical hydrostatic stress for PPSU, PSU, and PESU were determined to be able to compare the three systems together with literature results. The method used is based on the coupling of the evolution of yield stress and the evolution of embrittlement upon progressive aging as proposed by Engels et al. 5 For a specific annealing temperature, the measured time-to-embrittlement can be related to a value of the evolving yield stress, namely, a value of the state parameter S A . According to the method, for each annealing temperature, the value of the critical yield stress corresponding to embrittlement should, within experimental error, be the same. This was confirmed by Visser et al., 46 who showed for poly(vinyl chloride) (PVC) that the critical yield stresses corresponding to the ductile to brittle transition at three different anneal temperatures showed little variation. In other words, to a good approximation, there is a single critical thermodynamic state, reflected by a critical value of the state parameter, that marks the onset of embrittlement.
The evolution of embrittlement (top figures) and yield stress (bottom figures) upon aging for PSU and PESU is shown in Figure 5. The notch impact energy is calculated from the area under the force-displacement curve, divided by the crosssectional area behind the notch. The high energy levels correspond to ductile failure, manifested by large plastic deformation of the sample in the form of shear lips; the lower energy levels correspond to brittle failure, characterized by small to negligible plastic deformation and a rough fracture surface (see Figure 10). In full agreement with observations of Legrand 1 and Engels et al., 5 the transition from a ductile failure mode to a brittle one occurs on a shorter timescale for higher anneal temperatures. This transition is directly related to an increase in yield stress, shown in the lower figures. With increasing annealing time, the yield stress is observed to increase at all annealing temperatures. The vertical dashed lines in all figures represent the so-called time-to-embrittlement, taken as the time where the impact energy is midway between the ductile and brittle value, and are used to determine a corresponding critical yield stress (horizontal dashed lines). For PSU, these are 76.7, 77.4, and 76.7 MPa for anneal temperatures of 100, 120, and 140 C, respectively. For PESU, these were found to be 89.9, 90.5, and 90.7 for anneal temperatures of 100, 120, and 140 C respectively. These values are well within a range of 1 MPa, and in good agreement with the hypothesis of a single critical thermodynamic state dictating the transition from ductile to brittle failure. Additional evidence for the existence of a single critical thermodynamic state can be obtained by applying time-temperature superposition (TTS), which is shown in Figure 6 for PPSU. Using TTS, the data, both yield and embrittlement can be shifted to a similar reference temperature (in this case 150 C). Interestingly, both the evolution of yield stress as well as embrittlement can be shifted to a master curve using the same set of shift factors (depicted in Figure 7). An Arrhenius relation between shift factor and the reciprocal of the annealing temperature is obtained, which directly implies that a single activation energy applies to both the evolution of yield and embrittlement, in agreement with literature findings on polycarbonate. 5 The value of the critical yield stress is easily obtained and was found to be 77.6 MPa for PPSU.
The obtained critical yield stresses are a measure of the thermodynamic state of the material and are directly related to a unique value of the state parameter S A in the EGP model. 26 The critical state parameters are obtained with EGP model simulations of a tensile test, using a single axisymmetric element with uniaxial boundary conditions. The value of the state parameter is adapted to match the simulated yield stress to the experimentally determined critical yield stress, and hence the critical value S c a is obtained. This resulted in values of 24.9, 24.3, and 27.2 for PPSU-3, PESU, and PSU, respectively.
The obtained value of S c a is used as input for the EGP simulation of the impact experiment. In Figure 8, the result of the impact simulations for the three systems is shown. It must be noted that there is no failure criterion implemented in the model, implying that simulations are stopped after a certain amount of time to avoid excessive computation times. These times were chosen such that they were always well beyond  maximum hydrostatic stress occurring is extracted from each node on the symmetry line directly underneath the notch. Initially, the hydrostatic stress builds up close to the surface of the notch; however, with ongoing (plastic) deformation the maximum in hydrostatic stress follows the boundary of the plastic deformation zone toward the center of the bar. At the plastic zone boundary voiding is initiated, as was nicely demonstrated experimentally by Narisawa et al., 11 Gearing and Anand, 14 Figure 8 show images of the FE meshes at the point of reaching the critical hydrostatic stress. To show that the position of the maximum hydrostatic stress before failure corresponds to the position where voids are initiated, the simulation is compared to the fracture surface of a notch bar of PPSU-3 ( Figure 9). As can be seen, the positions indeed appear to match.
Interestingly, for PPSU-3, the maximum in hydrostatic stress is reached before a maximum in force. For PSU, the maximum in force is reached at the maximum occurring hydrostatic stress, whereas for PESU, the maximum in force is already reached before the maximum hydrostatic stress is reached. First, it should be noted that the hydrostatic stress reaches more of a plateau level rather than a sharp maximum. A possible explanation for the mismatch between the maximum in hydrostatic stress and experimentally observed maximum force for PPSU-3 is that the crazes that develop after voiding are still capable of carrying a load leading to a higher force before failure in experiments. It has been shown that this effect is more pronounced for high-molecular-weight materials, whereas for low-molecular-weight materials the maximum in hydrostatic stress and maximum in force appear to coincide. 5 For PESU, this obviously does not hold and a possible cause might be related to a rather different fracture surface than seen in PPSU and PSU, see Figure 10. For PESU, the cracks have propagated into a large volume behind the notch whereas PSU and PPSU only show crack propagation in the plane behind the notch. Another explanation could be the fact that a single critical yield stress, and thus state parameter is identified, whereas in reality a distribution of yield stress is observed over the crosssectional area of the tensile bar, especially for low mold temperatures. 47 Taking this distribution into account could yield improved simulation results. The critical value of the hydrostatic stresses obtained is 106.2 and 80.8 MPa for PPSU-3 and PSU respectively, indicated with the diamond markers in the figure. For PESU, the critical hydrostatic stress would be in that case be 90.4 MPa; however, because of the earlier occurrence of failure, the value is taken at the displacement-to-break observed in experiments and yields in that case 86.6 MPa.
In Figure 11, the obtained critical hydrostatic stresses are plotted as a function of entanglement density, including literature values for PC, 5 poly(methyl methacrylate) (PMMA), 12   PET, PPO, PMMA, and PSU are calculated from densities and molecular weights between entanglements reported by Fetters et al. 51 using: with ρ is the density, N a is Avogadro's constant, and M e is the molecular weight between entanglements. The entanglement density of PESU was directly derived from the plateau modulus measured with a dynamic mechanical thermal analysis (DMTA) experiment (ν e ¼ G 0 N =kT). The results in Figure 11 clearly show an increase in the critical hydrostatic stress with entanglement density which is in full agreement with literature. It can be concluded that the resistance to voiding increases with increasing entanglement density. It should be noted that all experiments and simulations in the current study have been performed at a temperature of 22 C, that is, the distance to the glass-transition temperature T g is different for each material. Comparing the hydrostatic stresses to T g does not yield a satisfactory trend. Deviations from the general trend could be caused by the fact that the critical hydrostatic stress depends on the molecular weight of the polymer: Engels et al., 5 but also Legrand,1 showed that molecular weight influences embrittlement. For a high-molecular-weight polymer, the transition from ductile to a brittle failure mode occurs on a later point in time than for a low-molecular-weight material, which directly implies an increase of the critical yield stress, and thus the critical hydrostatic stress, with increasing molecular weight. A similar observation was made by Pitman et al. 52 and Golden et al. 53 : a lower molecular weight yielded a lower craze stress. In contrast, the aging kinetics and intrinsic deformation response are hardly affected by molecular weight. 5,26 To show the effect of molecular weight, three other material grades of PPSU are examined in a similar way as PPSU-3. In Figure 12(a), the influence of   Figure 6. Solid line is a guide to the eye. molecular weight on embrittlement is shown. As can be seen, the time-to-embrittlement increases with increasing molecular weight. Furthermore, it appears that the initial value as well as the final value of the impact energy increase with increasing molecular weight. From a numerical point of view, since time-toembrittlement increases with molecular weight, the critical hydrostatic stress should increase likewise due to the increase of the critical state parameter. Numerical simulations are performed for each molecular weight with the value of the critical state parameter corresponding to the time-to-embrittlement of each grade. The results are summarized in Table II. The critical hydrostatic stresses found where 102.6, 103.8, and 105.8 MPa for PPSU-1, PPSU-2, and PPSU-4, respectively. The relation between number averaged molecular weight and critical hydrostatic stress is depicted in Figure 12(b). With an increase in molecular weight from ≈12 to ≈16 kg mol −1 , the critical hydrostatic stress increases approximately 3 MPa. The deviating value of PPSU-4, having the highest molecular weight but a similar critical stress as PPSU-3, could be caused by differences in polydispersity. Another reason could be the fact that M n is used, whereas Tervoort et al. 54 advocate the use of a corrected value: ϕM n * . Doing so, one corrects for the low-molecular-weight chains that actually dilute the system. Unfortunately, the molecular weight distributions were not available to perform this correction. Considering how well entangled all grades are, the ratio M w /M e ranges from 16 to 20, it is unlikely that a change from chain breakage to chain slip or vice versa occurs. It might very well be that the critical hydrostatic stress reaches a plateau value at high molecular weight, similar to what is often observed for the tensile strength 55 and fracture. 56

Predicting Embrittlement of Notched Tensile Bars
To further test the validity of a hydrostatic stress criterion, the aging induced embrittlement of notched tensile bars made from PPSU-3, having a different notch geometry (notch radius of 0.4 mm), is predicted using the obtained value of the critical hydrostatic stress for PPSU-3. The hypothesis is that the critical hydrostatic stress does not depend on notch geometry and that one can thus predict embrittlement using the critical hydrostatic stress as a criterion for the onset of brittle failure. EGP simulations are performed using again a mesh of an ASTM tensile bar, now with a notch radius of 0.4 mm. In Figure 13, the result for these simulations is shown for several values of the state parameter. Again, an increase in the state parameter yields an increase in the maximum hydrostatic stress. The gray dashed line represents the critical stress for PPSU-3 (σ c H ¼ 106.2 MPa). It was found that to reach the critical hydrostatic stress a value of 31.7 for the state parameter was needed; this value is clearly higher than obtained for a notch radius of 0.25 mm. Reason for this is the less severe build-up of hydrostatic stress behind the notch in case of a larger notch radius. This directly implies an increase in the time-to-embrittlement, that is, embrittlement occurs on a longer timescale. The value of the state parameter is used as input for a single element tensile simulation to determine the critical yield    Figure 14(b). For comparison, the data obtained for r = 0.25 from Figure 6(b) is added (light gray symbols). As expected, embrittlement of the notched tensile bars with a larger notch radius does indeed occur on a longer timescale. The ductile and brittle impact energy levels also seem to increase with increasing notch radius. The experimentally obtained time-to-embrittlement is indicated with the gray dashed line and was found to be 3.7 × 10 8 s, at a reference temperature of 150 C; in good agreement with the predicted value. The results confirm that the critical hydrostatic stress does not depend on notch geometry and proves to be a valid choice in predicting embrittlement.

CONCLUSIONS
Using a previously developed hybrid experimental-numerical method, the aging-induced embrittlement of three amorphous polymers was investigated. Coupling the evolution of yield stress of tensile bars to the evolution of embrittlement of notched tensile bars upon aging, the critical hydrostatic stress was determined for PPSU, PESU, and PSU. The obtained values correlate well with entanglement density: a denser network leads to a higher critical hydrostatic stress, that is, a higher resistance against voiding. Using the hydrostatic stress as a criterion the time-to-embrittlement of samples with a different notch geometry could be predicted. Although the critical hydrostatic stress appears to be material specific, it does depend on molecular weight. This given fact was investigated for PPSU and with an increase in molecular weight the critical hydrostatic stress, albeit in only a few megapascals, increases as well. All impact experiments and simulations in this study are performed at a single temperature and deformation rate. The influence of test temperature and deformation rate will be topic of a future publication.