A Study of Yielding and Plasticity of Rapid Prototyped ABS

Abstract

In this article, the yielding and plastic flow of a rapid-prototyped ABS compound was investigated for various plane stress states. The experimental procedures consisted of multiaxial tests performed on an Arcan device on specimens manufactured through photopolymerization. Numerical analyses were employed in order to determine the yield points for each stress state configuration. The results were used for the calibration of the Hosford yield criterion and flow potential. Numerical analyses performed on identical specimen models and test configurations yielded results that are in accordance with the experimental data.

1. Introduction

In recent years, rapid prototyping has experienced a significant growth, as the short product development times, corroborated with the fact that additive manufacturing equipment is becoming more accessible, has incentivized small enterprises to adopt such technologies for batch production [1,2].

Due to the different types of additive manufacturing technologies and parameters [2], the properties of the prototyped parts can show some variation [3,4]. For this reason, an extensive characterization is required in order to calibrate or develop material models that can be used in the numerical evaluation of the components during product design.

One of the main investigated topics of continuum mechanics from the 19th and 20th centuries was concerned with accurately describing the yielding behavior of materials for various spatial stress configurations and several phenomenological yield criteria were developed, based on experimental observations and yield hypotheses (Rankine—maximum principal stress, St. Venant—maximum principal strain, Tresca—maximum tangential stress, etc.) [5].

The most prevalent yield criterion used in engineering applications considers that yielding occurs when the distortional energy reaches a critical value. Though the exact origins are subject to debate (as J.C. Maxwell, T.M. Huber and H. Hencky addressed this hypothesis), the criterion is currently named after Richard von Mises [5].

Experimental results rarely if ever coincide with the theoretical yield surface predicted by the von Mises yield criterion. For some metals (such as steel alloys) the von Mises criterion determines better results than for others (such as some aluminum alloys) [6]. Other types of materials, such as soils and rocks, exhibit a completely different behavior and cannot be modeled with the von Mises criterion [7]. For this reason, a significant number of yielding models and flow potentials were developed throughout the years (Mohr–Coulomb, Drucker–Prager or Hill to name the most relevant).

Due to their wide range of compositions and microstructures, polymeric materials exhibit significant variations in terms of mechanical properties, and, in consequence, the corresponding yield criterion must be carefully chosen for each compound [8]. The von Mises yield criterion can determine accurate results for isotropic polymers, while the Hill criterion is more suitable for anisotropic compounds [9]. Some polymers exhibit different behavior when subjected to compression and tension, and thus, the Drucker–Prager model is recommended [10].

In consequence, the development and calibration of mathematical models is essential in accurately describing the mechanical behavior of polymeric materials. The numeric evaluation of mechanical properties is a mandatory part of product design, as complex geometries and loading conditions determine complex problems that cannot be solved analytically. Commercially available finite element analysis software implement a number of the above mentioned yield criteria, but they are usually limited to the most commonly used materials. As such, an important issue is the implementation of special material model routines in finite element codes, [11,12,13] especially for modeling plasticity and damage.

This study investigates the yielding behavior of a proprietary blend of ABS (acryonitrile–butadiene–styrene) specifically developed for photopolymerization applications. The proposed method of yield surface determination is based on tests performed on identical specimens for various stress configurations, using an Arcan device.

2. Experimental Procedures

The chosen geometry for the samples was based on the butterfly specimen, developed by Bai [6]. The specimen featured a circumscribed square with $50 \mathrm{mm}$ sides and two fillet radii ($R_1=33 \mathrm{mm}$ and $R_2=1.5 \mathrm{mm}$, Figure 1) in order to ensure stress concentration in the middle of the specimen and, consequently, a more facile identification of the first yield region.

The samples were manufactured through rapid prototyping using the PolyJet™ technology [14].

The Arcan device used in the experimental investigations allowed for seven orientations in increments of ${15}^{°}$. The pure shear loading was attributed for the orientation of $0^{°}$ and the plane strain tensile loading was attributed for the orientation of ${90}^{°}$ (Figure 2).

Preliminary numerical analyses using an isotropic linear elastic material model ($E=557 \mathrm{MPa}, \nu =0.3$) were performed in order to determine the relationship between the principal stresses for all the orientations. The geometric model (presented in Figure 1) was meshed using 43,634 s-order tetrahedral elements (C3D10) of varying sizes (denser mesh in the critical region). Beam-type multi-point constraints were generated between reference points (RPs) and the specimen clamping orifices. The specimen was rotated according to the testing angle while a displacement was imposed to one of the RPs, the other being fixed. The force-deflection curves were generated using the reaction from the fixed RP and the travel of the other RP.

The numerical results are plotted in Figure 3 in the principal stress space, along a von Mises yield surface for an easier visualization, and the relations between the principal stresses are presented in

Table 1. Relations between the principal stresses for the orientations of the Arcan device.

Orientation Angle	Principal Stress Relations	Orientation Angle	Principal Stress Relations
$0^{°}$	${\sigma }_1=-{\sigma }_2$	${60}^{°}$	${\sigma }_1=6.34{\sigma }_2$
${15}^{°}$	${\sigma }_1=-2.54{\sigma }_2$	${75}^{°}$	${\sigma }_1=4.88{\sigma }_2$
${30}^{°}$	${\sigma }_1=-10.46{\sigma }_2$	${90}^{°}$	${\sigma }_1=4.55{\sigma }_2$
${45}^{°}$	${\sigma }_1=14.98{\sigma }_2$

.

Considering that the results for the orientations at ${60}^{°}$, ${75}^{°}$ and ${90}^{°}$ yielded similar principal stress combinations, tests at ${60}^{°}$ and ${75}^{°}$ were omitted.

The resulting force-displacement curves are presented in Figure 4.

3. Yield Function Calibration

In order to determine the yield points for each orientation, numerical analyses were performed on the models described in the previous section and the resulting curves were compared to the experimental data. The material from the critical region of the specimen is considered to have yielded when the linear elastic response diverges from the experimental curve (Figure 5).

The crosshead travel corresponding to each yield point of a given specimen orientation was input as displacement for the RP in the numerical analyses, and the stress distribution was analyzed for the final step of the analysis, when yielding is considered to have occurred. Since the mid principal stress is negligible in the critical region (the surface from the middle of the specimen, where the two fillet radii meet, Figure 1), the volume elements that first yield can be considered to be in plane stress conditions. In consequence, the major and minor principal stresses were evaluated for the integration points corresponding to the critical region of the specimen (Figure 6), their combination determining the yielding point for the given orientation.

The resulting yield loci are plotted in Figure 7. It can be observed that the results do not coincide with the von Mises yield surface (Equation (1)).

$${\sigma }_y^2=3J_2$$

(1)

where ${\sigma }_y$ is the yield stress, $J_2$ is the second invariant of the deviatoric stress tensor ${\mathit{\sigma }}^{\prime }$,

$$J_2={\mathit{\sigma }}^{\prime }:{\mathit{\sigma }}^{\prime }=\frac{1}{2}{\displaystyle \sum }_{i,j=1}^3\left({\sigma }_{ii}^{\prime }{\sigma }_{jj}^{\prime }-{\sigma }_{ij}^{\prime }{\sigma }_{ji}^{\prime }\right)$$

(2a)

$${\mathit{\sigma }}^{\prime }=\mathit{\sigma }-{\sigma }_{Hyd}\mathit{I}$$

(2b)

$${\sigma }_{Hyd}=\mathrm{tr}\left(\mathit{\sigma }\right)$$

(2c)

$\mathit{\sigma }$ is the stress tensor and ${\sigma }_{Hyd}$ is the hydrostatic stress.

For a more accurate yield surface fit, the Hosford yield criterion was considered (Equation (3)) [15].

$$2{\sigma }_y^{2k}={\left({\sigma }_1-{\sigma }_2\right)}^{2k}+{\left({\sigma }_2-{\sigma }_3\right)}^{2k}+{\left({\sigma }_3-{\sigma }_1\right)}^{2k}$$

(3)

where $k$ is a material parameter and ${\sigma }_1$, ${\sigma }_2$ and ${\sigma }_3$ are the principal stresses (the eigenvalues of the stress tensor).

The comparison between the von Mises and the Hosford yield surfaces for $2k=6$ and $2k=8$ are presented in Figure 7 for the principal stress plane.

4. Multi-Linear Isotropic Hardening and Incremental Plasticity Model for the Hosford Criterion

The material was modeled as an isotropic elastic with multi-linear hardening. Isotropic elasticity was modeled using the generalized Hooke’s law (Equation (4)), with ${\mathit{ϵ}}^{\mathit{e}\mathit{l}}$ being the elastic strain tensor, $\mathit{a}$ the compliance tensor (fourth order symmetric tensor) and $\mathit{\sigma }$ the stress tensor.

$${\mathit{ϵ}}^{\mathit{e}\mathit{l}}=\mathit{a}:\mathit{\sigma }$$

(4)

For an isotropic material, in terms of the material constants, $E$ (the Young’s modulus) and $\nu$ (the Poisson ratio), the elastic strain tensor can be expressed as:

$${\mathit{ϵ}}^{\mathit{e}\mathit{l}}=\frac{1}{E}\left[\left(1+\nu \right)\mathit{\sigma }-\nu \mathit{I}\mathrm{tr}\left(\mathit{\sigma }\right)\right]$$

(5)

In a similar fashion, the stress tensor can be expressed as:

$$\mathit{\sigma }=2\mu \mathit{ϵ}+\kappa \mathit{I}\mathrm{tr}\left(\mathit{ϵ}\right)$$

(6)

where $\mu$ and $\kappa$ are the Lamé constants:

$$\mu =\frac{E}{2\left(1+\nu \right)}$$

(7a)

$$\kappa =\frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}$$

(7b)

Considering Equation (3), the yield function $f$ according to the Hosford criterion can be expressed as:

$$f={\left\{\frac{1}{2}\left[{\left({\sigma }_1-{\sigma }_2\right)}^{2k}+{\left({\sigma }_2-{\sigma }_3\right)}^{2k}+{\left({\sigma }_3-{\sigma }_1\right)}^{2k}\right]\right\}}^{\frac{1}{2k}}-{\sigma }_y\left({\overline{ϵ}}^{pl}\right)$$

(8)

where ${\overline{ϵ}}^{pl}$ is the equivalent plastic strain.

Plastic flow theory for associated models implies that there is no volume change during plastic deformation and thus, the components of the plastic strain rate tensor $\mathrm{d}{\mathit{ϵ}}^{pl}$ are proportional to the deviatoric stress tensor ${\mathit{\sigma }}^{\prime }$ through a scalar plastic multiplier $\mathrm{d}\lambda$:

$$\mathrm{d}{ϵ}_i^{pl}={\sigma }_i\mathrm{d}\lambda$$

(9)

which is known as the Levy–Mises flow rule [16].

The plastic volume consistency implies that the trace of the plastic strain rate tensor $\mathrm{d}{\mathit{ϵ}}^{pl}$ is null.

$$\mathrm{tr}\left(\mathrm{d}{\mathit{ϵ}}^{pl}\right)=\mathrm{d}{ϵ}_1^{pl}+\mathrm{d}{ϵ}_2^{pl}+\mathrm{d}{ϵ}_3^{pl}=0$$

(10)

Considering the relation between the total principal stress tensor components ${\sigma }_i$ and the principal deviatoric stress tensor components ${\sigma }_i^{\prime }$:

$${\sigma }_i={\sigma }_i^{\prime }+{\sigma }_{Hyd}$$

(11a)

$${\sigma }_{Hyd}=\frac{1}{3}\mathrm{tr}\left(\mathit{\sigma }\right)$$

(11b)

the Hosford equivalent stress can be expressed in terms of the principal deviatoric stress tensor components as:

$$\overline{\sigma }={\left\{\frac{1}{2}\left[{\left({\sigma }_1^{\prime }-{\sigma }_2^{\prime }\right)}^{2k}+{\left({\sigma }_2^{\prime }-{\sigma }_3^{\prime }\right)}^{2k}+{\left({\sigma }_3^{\prime }-{\sigma }_1^{\prime }\right)}^{2k}\right]\right\}}^{\frac{1}{2k}}$$

(12)

From Equation (9), one can assume that the equivalent plastic strain rate has a similar function to that of the equivalent stress for the Hosford yield criterion,

$$\mathrm{d}{\overline{ϵ}}^{pl}=C{\left\{\frac{1}{2}\left[{\left(\mathrm{d}{ϵ}_1^{pl}-\mathrm{d}{ϵ}_2^{pl}\right)}^{2k}+{\left(\mathrm{d}{ϵ}_2^{pl}-\mathrm{d}{ϵ}_3^{pl}\right)}^{2k}+{\left(\mathrm{d}{ϵ}_3^{pl}-\mathrm{d}{ϵ}_1^{pl}\right)}^{2k}\right]\right\}}^{\frac{1}{2k}}$$

(13)

where $\mathrm{d}{ϵ}_i^{pl}$ are the principal plastic strain rates (eigenvalues of the plastic strain rate ten-sor $\mathrm{d}{\mathit{ϵ}}^{pl}$) and $C$ is a constant that can be determined from the uniaxial loading case.

$$\left\{\begin{array}{c}\mathrm{d}{\overline{ϵ}}^{pl}=\mathrm{d}{ϵ}_1^{pl}\\ \mathrm{d}{ϵ}_2^{pl}=\mathrm{d}{ϵ}_3^{pl}=-\frac{1}{2}\mathrm{d}{ϵ}_1^{pl}\end{array}\right.$$

(14)

yielding $C=2/3$ and thus:

$$\mathrm{d}{\overline{ϵ}}^{pl}=\frac{2}{3}{\left\{\frac{1}{2}\left[{\left(\mathrm{d}{ϵ}_1^{pl}-\mathrm{d}{ϵ}_2^{pl}\right)}^{2k}+{\left(\mathrm{d}{ϵ}_2^{pl}-\mathrm{d}{ϵ}_3^{pl}\right)}^{2k}+{\left(\mathrm{d}{ϵ}_3^{pl}-\mathrm{d}{ϵ}_1^{pl}\right)}^{2k}\right]\right\}}^{\frac{1}{2k}}$$

(15)

The associated plastic flow hypothesis assumes that the yield surface expands through a plastic strain increment tensor that is normal to the tangent of the previous yield surface. In mathematical form, this is expressed as:

$$\mathrm{d}{\mathit{ϵ}}^{pl}=\mathrm{d}\lambda \frac{\partial f}{\partial \mathit{\sigma }}$$

(16)

The direction of the plastic strain increment is described by the tensor $\mathit{n}$,

$$\mathit{n}=\frac{\partial f}{\partial \mathit{\sigma }}$$

(17)

with the directions of the plastic strain increment being along the principal stress directions $n_i$:

$$\left\{\begin{array}{l}{n}_1=\frac{\partial f}{\partial {\sigma }_1}=\frac{1}{2}\left[{\left({\sigma }_1-{\sigma }_2\right)}^{\left(2k-1\right)}-{\left({\sigma }_3-{\sigma }_1\right)}^{\left(2k-1\right)}\right]{\overline{\sigma }}^{\left(1-2k\right)}\\ n_2=\frac{\partial f}{\partial {\sigma }_2}=\frac{1}{2}\left[{\left({\sigma }_2-{\sigma }_3\right)}^{\left(2k-1\right)}-{\left({\sigma }_1-{\sigma }_2\right)}^{\left(2k-1\right)}\right]{\overline{\sigma }}^{\left(1-2k\right)}\\ n_3=\frac{\partial f}{\partial {\sigma }_3}=\frac{1}{2}\left[{\left({\sigma }_3-{\sigma }_1\right)}^{\left(2k-1\right)}-{\left({\sigma }_2-{\sigma }_3\right)}^{\left(2k-1\right)}\right]{\overline{\sigma }}^{\left(1-2k\right)}\end{array}\right.$$

(18)

After some algebraic manipulations, it can be shown that:

$$\frac{\partial f}{\partial \mathit{\sigma }}=\frac{3}{2}\frac{{\mathit{\sigma }}^{\prime }}{\overline{\sigma }}$$

(19)

In incremental plasticity applications, it is worthwhile to link the scalar plastic multiplier $\mathrm{d}\lambda$ to the equivalent plastic strain $\mathrm{d}{\overline{ϵ}}^{pl}$ [16]. Considering the formulation of the Hosford yield criterion, the expression of $\mathrm{d}\lambda$ as a function of $\mathrm{d}{\overline{ϵ}}^{pl}$ is rather complex. Instead, numerical values were determined for several loading conditions such as uniaxial (${\sigma }_1={\sigma }_y, {\sigma }_2={\sigma }_3=0$), equibiaxial (${\sigma }_1={\sigma }_2, {\sigma }_3=0$), proportional biaxial (${\sigma }_1=\alpha {\sigma }_2, {\sigma }_3=0$ with $\alpha$ being an integer value), pure shear (${\sigma }_1=-{\sigma }_2, {\sigma }_3=0$), etc. It was observed that $\mathrm{d}\lambda =\left(0.93..1.07\right)\mathrm{d}{\overline{ϵ}}^{pl}$ and thus it was considered that:

$$\mathrm{d}\lambda \cong \mathrm{d}{\overline{ϵ}}^{pl}$$

(20a)

$$\mathrm{d}{\mathit{ϵ}}^{pl}=\mathrm{d}{\overline{ϵ}}^{pl}\frac{3}{2}\frac{{\mathit{\sigma }}^{\prime }}{\overline{\sigma }}$$

(20b)

Based on the incremental form of elasticity and the strain decomposition principle:

$${\mathit{ϵ}}^{tot}={\mathit{ϵ}}^{el}+{\mathit{ϵ}}^{pl}$$

(21)

where ${\mathit{ϵ}}^{tot}$ is the total strain tensor, ${\mathit{ϵ}}^{el}$ is the elastic strain tensor and ${\mathit{ϵ}}^{pl}$ is the plastic strain tensor, the elastic strain tensor can be expressed as:

$${\mathit{ϵ}}_{t+1}^{el}={\mathit{ϵ}}_t^{el}+\Delta {\mathit{ϵ}}^{el}={\mathit{ϵ}}_t^{el}+\Delta {\mathit{ϵ}}^{tot}-\Delta {\mathit{ϵ}}^{pl}$$

(22)

and the stress tensor can be expressed as:

$$\mathit{\sigma }=2\mu \left({\mathit{ϵ}}_t^{el}+\Delta {\mathit{ϵ}}^{tot}\right)+\kappa ·\mathrm{tr}\left({\mathit{ϵ}}_t^{el}+\Delta {\mathit{ϵ}}^{tot}\right)\mathit{I}-2\mu \Delta {\mathit{ϵ}}^{pl}$$

(23)

The first part of Equation (23) is called the elastic predictor stress,

$${\mathit{\sigma }}^{pr}=2\mu \left({\mathit{ϵ}}_t^{el}+\Delta {\mathit{ϵ}}^{tot}\right)+\kappa ·\mathrm{tr}\left({\mathit{ϵ}}_t^{el}+\Delta {\mathit{ϵ}}^{tot}\right)\mathit{I}$$

(24)

and the second part represents a plastic corrector.

Considering Equation (20b), the stress tensor can be expressed as:

$$\mathit{\sigma }={\mathit{\sigma }}^{pr}-2\mu \Delta {\overline{ϵ}}^{pl}\frac{3}{2}\frac{{\mathit{\sigma }}^{\prime }}{\overline{\sigma }}$$

(25)

From Equation (19b), Equation (25) can be written in terms of the deviatoric stress tensor as:

$${\mathit{\sigma }}^{\prime }\left(1+3\mu \frac{\Delta {\overline{ϵ}}^{pl}}{\overline{\sigma }}\right)={\mathit{\sigma }}^{pr}-\frac{1}{3}\mathrm{tr}\left(\mathit{\sigma }\right)\mathit{I}$$

(26)

with the second term of Equation (26) being the deviatoric stress predictor:

$${\mathit{\sigma }}^{pr}-\frac{1}{3}\mathrm{tr}\left(\mathit{\sigma }\right)\mathit{I}={\left({\mathit{\sigma }}^{\prime }\right)}^{pr}$$

(27)

In order to evaluate if the material has reached the yield point, an equivalent predictor stress ${\overline{\sigma }}^{pr}$ is introduced:

$${\overline{\sigma }}^{pr}={\left\{\frac{1}{2}{\left({\sigma }_1^{pr}-{\sigma }_2^{pr}\right)}^{2k}+{\left({\sigma }_2^{pr}-{\sigma }_3^{pr}\right)}^{2k}+{\left({\sigma }_3^{pr}-{\sigma }_1^{pr}\right)}^{2k}\right\}}^{\frac{1}{2k}}$$

(28)

In consequence, Equation (26) must also be expressed using equivalent values:

$$\overline{\sigma }\left(1+3\mu \frac{\Delta {\overline{ϵ}}^{pl}}{\overline{\sigma }}\right)={\overline{\sigma }}^{pr}$$

(29a)

$$\overline{\sigma }={\overline{\sigma }}^{pr}-3\mu \Delta {\overline{ϵ}}^{pl}$$

(29b)

For the plastic flow law, an isotropic multi-linear hardening model was considered. The multi-linear hardening model assumes that the equivalent plastic stress is a discrete function of the equivalent plastic strain, the behavior being defined with the help of linear interpolations between the input points (Figure 8).

Since no uniaxial tests were performed on this material, the hardening curve was approximated from the tests performed at ${45}^{°}$ [17] and fine-tuned using finite element analyses.

The isotropic hardening model assumes that the yield surface expands with the accumulation of plastic strain through the help of a hardening function $H\left({\overline{ϵ}}^{pl}\right)$.

$$f={\overline{\sigma }}^{pr}-3\mu \Delta {\overline{ϵ}}^{pl}-H\left({\overline{ϵ}}^{pl}\right)-{\sigma }_y=0$$

(30)

Relation (30) is a non-linear equation with respect to the equivalent plastic strain increment $\Delta {\overline{ϵ}}^{pl}$. The implicit integration method considers the Newton–Raphson algorithm in solving the equation [18]. Expanding the yield function $f$ in a Taylor series and truncating after the first term yields:

$$f+\frac{\partial f}{\partial \Delta {\overline{ϵ}}^{pl}}d\left(\Delta {\overline{ϵ}}^{pl}\right)=0$$

(31)

with

$$\frac{\partial f}{\partial \Delta {\overline{ϵ}}^{pl}}=-3\mu -\frac{\partial H\left({\overline{ϵ}}^{pl}\right)}{\partial \Delta {\overline{ϵ}}^{pl}}$$

(32)

The linear hardening function can be expressed in rate form as:

$$H\left({\overline{ϵ}}^{pl}\right)=h\left({\overline{ϵ}}^{pl}\right)·\Delta {\overline{ϵ}}^{pl}$$

(33)

where $h\left({\overline{ϵ}}^{pl}\right)$ is the slope of the hardening curve (the plastic modulus):

$$h\left({\overline{ϵ}}^{pl}\right)=\frac{\Delta {\sigma }^y}{\Delta {\overline{ϵ}}^{pl}}$$

(34)

For the multi-linear hardening model, the hardening parameter $h$ is evaluated from the yield stress-plastic strain values (${\sigma }_i^y,{\overline{ϵ}}_i^{pl}$):

$$h_i=\frac{{\sigma }_{i+1}^y-{\sigma }_i^y}{{\overline{ϵ}}_{i+1}^{pl}-{\overline{ϵ}}_i^{pl}}$$

(35)

In consequence, the derivative of the equivalent plastic strain increment $d\left(\Delta {\overline{ϵ}}^{pl}\right)$ is obtained by introducing Equations (30), (32) and (33) in Equation (31):

$$d\left(\Delta {\overline{ϵ}}^{pl}\right)=\frac{{\overline{\sigma }}^{pr}-3\mu \Delta {\overline{ϵ}}^{pl}-H\left({\overline{ϵ}}^{pl}\right)-{\sigma }^y\left(0\right)}{3\mu +h_i}$$

(36)

The hardening function can be expressed in incremental form as:

$$H_{t+1}=H_t+h_t\Delta {\overline{ϵ}}_t^{pl}$$

(37)

and the plastic strain increment as:

$${\overline{ϵ}}_{t+1}^{pl}=\Delta {\overline{ϵ}}_t^{pl}+d\left(\Delta {\overline{ϵ}}^{pl}\right)$$

(38)

5. Numerical Results and Discussion

The implicit integration scheme presented in Paragraph 5 was implemented in the commercial software Abaqus 6.14 using a UMAT subroutine [19]. The numerical analyses were re-performed using a von Mises material and the Hosford model for $2k=8$, both formulations using the hardening curve described in the previous paragraph. The results are presented in Figure 9 for all the configurations of the Arcan device.

As a general observation, apart from the fact that the Hosford criterion determines lower yield point values (as predicted, Figure 7), the von Mises flow potential determines higher stress values for similar equivalent plastic strains. This aspect is more pronounced at the rotation angles of $0^{°}$ and ${15}^{°}$, which are associated with the stress combinations for which the two yield functions show the biggest discrepancies (the elements are predominantly subjected to shear loading).

In addition, it can be observed that, for the stress states that lead to large plastic deformations before failure (Arcan device orientations at $0^{°}$ and ${15}^{°}$, Figure 9a,b, respectively), both models overestimate the stress response after a certain equivalent plastic strain. This aspect may be caused by damage onset in the material, which can progressively reduce the stiffness of the material [20]. Figure 10 shows an image capture from the tests performed at ${15}^{°}$ orientation, where it was observed that, prior to the specimen failure, a void nucleation can be observed in the critical region, which reduces the effective cross-section area and eventually causes a slight drop in the reaction force in the plateau region (Figure 9b).

The failure mechanism associated with ductile materials considers that the critical plastic strain at damage onset is a function of the stress triaxiality $\eta$, the Lode angle parameter $\mathsf{\Theta }$ and the equivalent plastic strain rate ${\dot{\overline{ϵ}}}^{pl}$ [20].

$$\eta =\frac{{\sigma }_{Hyd}}{\sqrt{3J_2}}$$

(39a)

$$\mathsf{\Theta }=\arccos \left(\frac{3}{2}\sqrt{3}\frac{J_3}{{J_2}^{\frac{3}{2}}}\right)$$

(39b)

$$J_3=\det \left({\mathit{\sigma }}^{\prime }\right)$$

(39c)

where $J_3$ is the third invariant of the deviatoric stress tensor. This polymer compound exhibits higher critical plastic strains at lower stress triaxiality values (associated with the test orientations at $0^{°}$, ${15}^{°}$ and ${30}^{°}$), and the inaccuracies in modeling of high plastic strains may be a shortcoming of the associated flow rule that was considered.

6. Conclusions

This work proposes the implementation of the Hosford criterion in modeling the plasticity of rapid prototyped polymers. Experimental procedures were performed in order to determine the yield surface of the investigated polymer. Numerical analyses were corroborated with the experimental data to calibrate the Hosford model, and the yield function and flow potential were implemented in the commercial software Abaqus. The Hosford criterion and flow potential determined better results than the von Mises model, especially for the tests where shear had a significant effect (orientations of $0^{°}$, ${15}^{°}$ and ${30}^{°}$). The implemented model proved very accurate for relatively low equivalent plastic strains ($<0.5 \mathrm{mm}/\mathrm{mm}$), higher plastic strains determining an overestimation in reaction force for the numerical model (peak values of 6.4% for the tests at $0^{°}$ and 5.5% for the tests at ${15}^{°}$).

Overall, the results prove satisfactory, but future work will focus on improving the flow potential in order to increase the numerical analyses’ accuracy. A non-associative flow rule will be considered in order to obtain better results at high plastic strains.

Considering that the materials used in rapid prototyping applications are usually provided by the equipment manufacturer and consist of proprietary blends, there is no control over the chemical properties and little control over the mechanical properties (through the tweaking of manufacturing parameters). In consequence, each material compound must be investigated individually and modeled according to its mechanical properties in order to ensure accurate product analysis results.