Linear Drucker-Prager Yield Criterion Calibration for Polypropylene under Multiaxial Load

1. Introduction

The introduction of polymeric materials, specifically polypropylene, to the market in 1951 [1], represented a pivotal advancement in the field of materials science and engineering. These polymers have since become integral across a diverse array of industrial applications like automotive, medical, construction, electronics and consumer goods [2,3,4]. The preference for polymeric materials, particularly polypropylene, is based on their outstanding mechanical properties, intrinsic resistance to aging, capability to maintain dimensional stability throughout their lifecycle, superior resistance to crack initiation or propagation, and their suitability for high-volume manufacturing processes such as extrusion and injection molding [5,6,7].

Effective design process of designing and validating component made of polymers in industries requires the usage of simulation methods, including the Finite Element Method [8]. To make this possible, it is necessary to accurately represent the behavior of plastics under the different loading regimes and to relative them in the numerical material model [9,10,11]. In addition to the commonly known nonlinear behavior of thermoplastics, influence of strain rate, temperature, stress state and loading type [12,13,14,15], the critical phenomena is the influence of hydrostatic pressure on the yield stress, plastic flow, the destruction criterion [16,17,18]. Studies conducted in the past have clearly demonstrated the occurrence of different elastoplastic properties for loads of different characteristics resulting from tests ranging from compression, through shear, to tension [19,20,21,22].

One of the widely used constitutive model capable to reproduce the plastic behavior is the linear Drucker-Prager criterion, which makes the transition of the material to the plastic range, i.e. the value of the yield stress, dependent on the value of hydrostatic pressure. Thanks to this, it is possible to distinguish the behavior of polymers under different load conditions, including fact that for most of thermoplastics the yield strength for compression is higher than for tension [23,24,25,26]. The method originally comes from soil research but it has found wide application for plastics, composites and metals [27,28,29,30].

This article focuses on the use of the Arcana machine which is widely used when the only possible form of material samples are flat plates of a defined thickness [31,32,33,34,35,36]. It is necessary to properly cut the plates to obtain a multiaxial state of tension using an appropriate fillet radius at the necking point. An additional verification is a uniaxial compression test with extend criteria to compression state. The obtained material model is finally tested in the Abaqus numerical software to confirm the compliance of the obtained parameters for the Drucker-Prager model and the numerical simulations carried out using it with the actual behavior of polypropylene samples.

2. Material Behavior

Stress tensor (${\mathit{\sigma }}_{\mathit{i}\mathit{j}}$) according to the matrix form from Equation 1 could be explain using hydrostatic (mean) (${\mathit{\sigma }}_{\mathit{h}}$) stress tensor responsible for volumetric change and deviatoric stress tensor (${\mathit{s}}_{\mathit{i}\mathit{j}}$) cause body shape deviation.

$${\sigma }_{ij}=\left(\begin{array}{ccc}{\sigma }_{xx}& {\tau }_{xy}& {\tau }_{xz}\\ {\tau }_{yx}& {\sigma }_{yy}& {\tau }_{yz}\\ {\tau }_{zx}& {\tau }_{zy}& {\sigma }_{zz}\end{array}\right)=\left(\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{h}}& 0& 0\\ 0& {\mathit{\sigma }}_{\mathit{h}}& 0\\ 0& 0& {\mathit{\sigma }}_{\mathit{h}}\end{array}\right)+\left(\begin{array}{ccc}{s}_{xx}& s_{xy}& s_{xz}\\ s_{yx}& s_{yy}& s_{yz}\\ s_{zx}& s_{zy}& s_{zz}\end{array}\right)={\sigma }_h{\delta }_{ij}+s_{ij}$$

(1)

where Kronecker Delta function is defined by:

$${\delta }_{ij}=\left\{\begin{array}{c}1,ifi=j\\ 0,ifi\ne j\end{array}\right.$$

(2)

Dilatation state defined by pressure ($\mathit{p}$) could be drawn using stress tensor first invariant (${\mathit{I}}_1$):

$$p=-{\sigma }_h=-{\displaystyle \frac{1}{3}}trace\left({\sigma }_{ij}\right)=-{\displaystyle \frac{{\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz}}{3}}=-{\displaystyle \frac{I_1}{3}}$$

(3)

Equivalent stress called von Mises ($\mathit{q}$) is defined by second invariant of deviator state (${\mathit{J}}_2$):

$$q={\sigma }_{eq}=\sqrt{3{\mathit{J}}_2}$$

(4)

$${\mathit{J}}_2={\displaystyle \frac{{\left({\sigma }_{xx}-{\sigma }_{yy}\right)}^2+{\left({\sigma }_{yy}-{\sigma }_{zz}\right)}^2{+\left({\sigma }_{zz}-{\sigma }_{xx}\right)}^2}{6}}+{\tau }_{xy}+{\tau }_{yz}+{\tau }_{zx}$$

(5)

Thirst invariant of deviator state (${\mathit{J}}_3$) needed to explain Drucker-Prager criterion could by defined by:

$${r=J}_3=det\left(s_{ij}\right)$$

(6)

$$det\left(s_{ij}\right)={\sigma }_{xx}{\sigma }_{yy}{\sigma }_{zz}+{\tau }_{xy}{\tau }_{yz}{\tau }_{zx}+{\tau }_{xz}{\tau }_{yx}{\tau }_{zy}-{\tau }_{xz}{\sigma }_{yy}{\tau }_{zx}-{\tau }_{xy}{\tau }_{yx}{\sigma }_{zz}-{\sigma }_{xx}{\tau }_{yz}{\tau }_{zy}$$

(7)

The linear Drucker–Prager criterion commonly used for polymers can be represented in the meridian plane by the yield surface using deviatoric stress factor ($\mathit{t}$) as a function of hydrostatic pressure ($\mathit{p}$) as in Figure 1, where the parameters defining its location are the friction angle ($\mathit{\beta }$) and material cohesion ($\mathit{d}$) according to Equation 8:

$$t=d+\mathrm{t}\mathrm{a}\mathrm{n}\left(\beta \right)$$

(8)

where deviatoric stress factor ($\mathit{t}$) could be defined by Equation 9:

$$t={\displaystyle \frac{1}{2}}q\left(1+{\displaystyle \frac{1}{K}}-\left(1-{\displaystyle \frac{1}{K}}\right){\left({\displaystyle \frac{r}{q}}\right)}^3\right)$$

(9)

using ratio of yield stress for triaxial tension and compression ($\mathit{K}$), what could be described on deviatoric plane in Figure 2. Based on the research on Polypropylene [22] and due to the limited scope of possible tests on the Arcana machine, a K factor of 1 was assumed. Due to this $\mathit{t}=\mathit{q}$ on Figure 1.

Plastic flow direction ($\mathit{d}{\mathit{\epsilon }}^{\mathit{p}\mathit{l}}$) for polymer could be defined in Equation 10 by dilatation angle ($\mathit{\psi }$) using Poison ration ($\mathit{\nu }$) [37,38].

$$\psi ={tan}^{-1}\left({\displaystyle \frac{3\left(1-2\nu \right)}{2\left(1+\nu \right)}}\right)$$

(10)

In Equation 8:

• tangent of friction angle ($\mathit{\beta }$) represent the sensitivity for hydrostatic pressure ($\mathit{p}$),
• material cohesion ($\mathit{d}$) reflect isotropic hardening, which could be defined by Equation 11 using uniaxial tension (${\mathit{\sigma }}_{\mathit{t}}$), compression (${\mathit{\sigma }}_{\mathit{c}}$) or shear ($\mathit{\tau }$) yield stress,

  $$d=\left({\displaystyle \frac{1}{K}}+{\displaystyle \frac{1}{3}}tan\left(\beta \right)\right){\mathit{\sigma }}_t=\left(1-{\displaystyle \frac{1}{3}}tan\left(\beta \right)\right){\mathit{\sigma }}_{\mathit{c}}={\displaystyle \frac{\sqrt{3}}{2}}\left(1+{\displaystyle \frac{1}{K}}\right)\mathit{\tau }$$

  (11)

Drucker-Prager yield surface in the reverse principal axes is a cone around hydrostatic axis with a circular cross section for K factor equal to 1 as on Figure 3.

To clearly determine the loading state of the material, the triaxiality ($\mathit{\eta }$) parameter (Equation 12) should be used, the basic values of which are given in Table 1.

$$\eta =-{\displaystyle \frac{p}{q}}$$

(12)

3. Material, Samples and Experimental Techniques

The following tests were performed for thermoplastic Sabic 83MF10 material, universal polypropylene polymer provided by Saudi Basic Industries Corporation (SABIC) and the basic parameters of which are presented in Table 2. Dilatation angle was calculated based on Equation 2.

Due to the fact that despite the lack of reinforcing inclusions in the material such as fibers, noticeable anisotropy of properties is possible [41,42], it was decided to cut samples from the plate far away from the injection points. This allowed to obtain a quasi-random arrangement of molecules.

Tension and shear samples subjected to analysis on a Zwick type universal testing machine with Arcana handle was obtained in the form of flat plates received in injection molding process, with a thickness of $\mathrm{2,2}\left[\mathit{m}\mathit{m}\right]$. The samples were machine cut with the radius of $1.5\left[\mathit{m}\mathit{m}\right]$ and dimensions described on Figure 4.

The sample for uniaxial compression tests is a single cylinder with a diameter of $3\left[\mathit{m}\mathit{m}\right]$ cut from the previously mentioned plates.

To achieve wide range of triaxiality $\left(\mathit{\eta }\right)$ and according to it different load type in stated between tension and shear, the Arcana machine was used (Figure 5). This type of holder allows to stretch samples at angles from 0 to 90 degrees. It consists of a pair of holders that secure it to the testing machine, two semicircular discs with holes that allow for mounting at 7.5 degrees and holders that grip the tested sample directly. The entire structure is made of aluminum alloy 7075 (Fortal), ensuring high rigidity of the structure.

All test on Arcana machine were performed with constant head travel speed of $1\left[\mathit{m}\mathit{m}/\mathit{s}\right]$ which allows to achieve strain rate ($\dot{\mathit{\epsilon }}$) between $0.06$ and $0.11\left[1/\mathit{s}\right]$ during the tests. Due to the selected machine speed it was possible to avoid, on the one hand, the influence of material hardening at high strain rates described e.g. by Cowper-Symonds model [39,40] and, on the other hand, the effect of polymer relaxation [43,44].

4. Experimental Data Processing and Results

The linear Drucker-Prager criterion was chosen for further calibration due to the impossibility of obtaining results for triaxial tension on flat samples with undercut. Figure 6 shows the results of tension (and shear) tests on the Arcana machine and the result of the compression test of a cylindrical sample. For each variant, three tests were performed, and two extreme cases (max and min values) in terms of yield strength were rejected from further evaluation.

The following tests were performed to determine the yield strength, which according to ISO 527-1 defines yielding as the first deviation in true stress-strain curve from pure elastic response due to the lack of a distinct yield strength. Force ($\mathit{F}$) and displacement ($\mathit{d}\mathit{i}\mathit{s}\mathit{p}$) were measure using Zwick embedded cross-head sensor. The values of triaxiality ($\mathit{\eta }$), pressure ($\mathit{p}$) and equivalent von Mises stress ($\mathit{q}$) were determined from FEM simulations representing the behavior of the samples based on the calibration of the material curve with the load case as a average value for whole cross section in the critical most strained surface. All results needed for further evaluation were described below in Table 3.

The larger the angle on the Arcana machine, the lower the yield strength seems to be. This does not reflect the actual parameters, due to the complex stress state, different for each sample. The decreasing apparent stiffness of the sample with increasing angle results from the transition from pure tension to pure shear, where instead of Young's modulus ($\mathit{E}$), the basic parameter defining elasticity is the Kirchhoff’s ($\mathit{G}$) modulus, which calculated according to Equation 13 is 2.8 times lower.

$$G={\displaystyle \frac{E}{2\left(1-\nu \right)}}$$

(13)

All the result in x relationship were described as a linear function of Drucker-Prager yield criteria in Figure 7. The regression line was obtained using solver which has the objective function with scope to minimize the ${\mathit{R}}^2$ error value.

As can be seen in the above graph, the results obtained for the Arcana machine correspond to the values from pure shear ($\mathit{\eta }=0$) to biaxial tension ($\mathit{\eta }=2/3$), which was possible thanks to the proposed undercuts on flat samples. The disc compressed on the strength machine obtained a result slightly lower ($\mathit{\eta }<-1/3$) than for uniaxial compression due to its relatively small thickness and the friction occurring in the test between the upper and lower external surfaces of the sample and the contact linings of the testing machine. Final values for linear Drucker-Prager criteria were described in Table 4.

All points achieved an average ${\mathit{R}}^2$ error value of $0.35$, which suggests a very good fit to the linear relationship. An extension of the proposed linear criterion can be the tension cut-off, which determines the limiting value of the tensile pressure (${\mathit{p}}_{\mathit{t}}$), beyond which plastic deformation occurs, regardless of the value of von Mises stress ($\mathit{q}$) [45].

In order to confirm the correct validation of the Drucker-Prager model, simulations were performed on previously tested samples and the calculation results in Abaqus were compared to the values obtained on the strength test bench. Implicit simulations were used using Abaqus Standard solver. The discretized models used in the simulations (shown in Figure 8) for Arcana and compression have respectively 16263 C3D10, 10-node quadratic tetrahedron elements with characteristic size of $0.5\left[\mathit{m}\mathit{m}\right]$ in critical space and 2200 C3D20, 20-node quadratic hexahedral elements with $0.2\left[\mathit{m}\mathit{m}\right]$ for whole geometry.

The finite element simulation was performed in a way that was designed to reproduce the actual test conditions as precisely as possible. For the samples tested on the Arcana machine, the model was completely fixed on the lower side in the locations indicated as mounting place in Figure 4. On the other side, a displacement was imposed in the vertical direction. The Arcana machine itself was not included in the simulation due to the stiffness of the entire structure made of aluminum that was several orders of magnitude higher. For the disc compressed on both sides, contact surfaces were simulated with a friction coefficient of $0.2\left[-\right]$, of which only one, the upper surface, moves towards the lower fixed plane. In order to follow the simulation progress above the yield point, the yield curve shown in Figure 9 was determined using numerical calibration from sample no 1 (Arcana – 0deg) test curve.

The hardening behavior was determined for range from 0 to 0.4 of plastic strain and implemented in Abaqus numerical software in tabular form with 17 points. Based on the above data, a comparison (on Figure 10) between previous laboratory results and the numerical simulation results was made for all proposed samples.

As shown in the above graph, the simulation result shows the moment of yielding with a satisfactory result. The moment of appearance of the first plasticized finite elements and the transition of the entire critical cross-section into the plastic range are clearly visible. Due to the determination of the plasticity curve used for simulation based on sample no. 1, the simulation course outside the elastic range is best shown. For models with an angle of setting above 30°, the laboratory results are characterized by a less clear moment of yielding compared to the simulation results. This may result from the settings in the Abaqus program, which do not take into account the possibility of small angular movements of the samples at the place of attachment, which can occur during tests on the Arcana machine. For the compressed sample, the moment of yielding is shown satisfactorily in the numerical model, however, a comparison of the course in the plastic range shows clearly that the yield curve differs significantly from that obtained from the tensile sample. The anisotropic nature of the elastic range is noticeable, which is demonstrated by the fact that for Arcana machine settings close to 90° introducing a shear state in the sample, the Young Modulus ($\mathit{E}$) values ​​are nearly 10% lower than those originally used in the simulation.

5. Conclusions and Areas for Further Research

The comprehensive study of Drucker-Prager presented in this article has shown a noticeable and significant linear dependence of the yield strength on the material state and, consequently, an unequivocal effect of pressure on the criterion under study. The proposed approach using a specially manufactured Arcana machine with samples with special undercuts allows for the representation of the state from pure shear to biaxial tension, which is not possible using standard samples. The supplementation and extension of the studies to the scope of uniaxial compression is a compression test performed on a slender disc. The obtained results uniformly achieve a linear character of variation.

All previously laboratory tested samples were numerically simulated, which allowed us to confirm the ability of the Drucker-Prager criterion implemented in the Abaqus program to properly reflect the moment of plasticization of samples. Due to the demonstrated variable course of the plasticization curve for different material loading states (Triaxiality), it seems advisable to take this phenomenon into account in future studies. At the moment, there are no ready-made solutions. Due to the stiffening of the numerical model of the sample mounting points in the holder, a small deviation occurred for the elastic range in samples with a machine inclination angle >15°.

Further research should focus on the longer range of Triaxiality, above 2-axis stretching and determining the so-called tension cut-off, which defines the pressure at which the material immediately plasticizes. A certain simplification of this phenomenon may be an attempt to calibrate the exponent or hyperbolic model, which indirectly takes into account the mentioned phenomenon.