*2.2. Graded Tensile Tests*

Graded tensile tests were carried out at three different temperatures (44 ◦C, 60 ◦C, and 80 ◦C), with the same deformation step size of 0.25 mm (see Figure 3e). The specimen was elongated by 0.25 mm at each step, at a deformation rate of 0.017 mm s<sup>−</sup>1. The time delay at the given strain value was always 60 s. This was done until the specimen failed. During the graded tensile test, the creep phenomenon was partially detected. Up to the value of the tensile force of approximately 800 N, there was no significant relaxation of the tension.

### *2.3. Indentation Tests*

Indentation tests were carried out at 23 ◦C, 44 ◦C, 60 ◦C, and 80 ◦C (see Figure 3f) on the samples depicted in Figure 2. The area where the indentation test was carried out is shown by a dashed rectangle. The indenter was a steel sphere with a diameter of 5 mm. The indentation test consisted of three phases. First, the indenter was pressed into the material at a speed of 0.017 mm s−<sup>1</sup> to a depth of 0.5 mm, then the time delay of 300 s followed, and finally the indenter returned to the starting position (relief) at the same speed.

#### **3. Material Model**

One of the material models used for viscoplastic materials is the Anand material model [12]. The Anand model was proposed for use in the analysis of the rate-dependent deformation of metals at high temperatures. The Anand viscoplastic model is pre-built in the commercial finite element software ANSYS©, and therefore it is much easier to use. The Anand model is typically used for solder alloys [21,22], however the previously mentioned benefits led the authors to test it for materials used for 3D printing (ABS-M30). The material model involves 11 material parameters: the Poisson ratio (*u*), Young's modulus (E), initial value of deformation resistance (s0), activation energy/universal gas constant (Q/R), pre-exponential factor (A), stress multiplier (xi), strain rate sensitivity of stress (m), hardening/softening constant (h0), coefficient for deformation resistance saturation value (S), strain rate sensitivity of saturation (deformation resistance) value (n), and strain rate sensitivity of hardening or softening (a).

The Anand viscoplastic model was originally developed for material forming applications [12,13]. It is also applicable to general viscosity problems, which include the influence of the strain rate and temperature. Materials at elevated temperatures are highly dependent on the influence of the temperature magnitude and history, strain rate, and strain hardening. The Anand model is a complex material model that introduced an internal variable S (deformation resistance), a variable that represents the resistance against the plastic behavior of the material. This is different from other material models, and the decomposition of the individual components of the deformation is not straightforward.

The rate of plastic deformation is described using the following relationship:

$$
\dot{\varepsilon}\_{\rm pl} = \dot{\varepsilon}\_{\rm pl}^{\rm a} \left( \frac{3}{2} \frac{\rm S}{\rm q} \right),
\tag{1}
$$

where **.** <sup>ε</sup>pl is the tensor of the inelastic strain rate and . ε a pl is the rate of accumulated equivalent plastic strain. Here, . ε a pl is given by the equation:

$$
\dot{\varepsilon}^{\text{a}}\_{\text{pl}} = \left(\frac{2}{3} \dot{\varepsilon}\_{\text{pl}} : \dot{\varepsilon}\_{\text{pl}}\right)^{\frac{1}{2}},\tag{2}
$$

where the operator ":" stands for inner product of the tensors. S is the deviator of the Cauchy stress tensor, which can be expressed using the following relation:

$$\mathbf{S} = \mathbf{\sigma} - \text{pl},\tag{3}$$

where σ is the Cauchy stress tensor and p is defined as one-third of the trace of the tensor matrix σ, as in the following relation:

$$\mathbf{p} = \frac{1}{3}\text{tr}(\boldsymbol{\sigma}).\tag{4}$$

Here, I represents a second-order unit tensor. The quantity q is the equivalent stress according to the following relation:

$$\mathbf{q} = \left(\frac{\mathbf{3}}{2}\mathbf{S} : \mathbf{S}\right)^{\frac{1}{2}}.\tag{5}$$

The rate of accumulated plastic deformation depends on q and on the internal state variable s. This dependence can be expressed by the following equation:

$$\dot{\varepsilon}^{\text{a}}{}\_{\text{pl}} = \text{A} \mathbf{e}^{(-\frac{\mathbf{Q}}{\text{BT}})} \left\{ \sinh \xi \frac{\mathbf{q}}{\text{s}} \right\}^{\frac{1}{m}} \text{.} \tag{6}$$

where A, ξ, and m are the model constants; Q is the activation energy; R is the universal gas constant; θ is the absolute temperature; and s is the internal state variable.

Equation (6) indicates that plastic deformation occurs at any stress level. This contrasts with other theories of plasticity that use areas of plasticity (yield function). Classical theories assume that plasticity occurs above a certain stress value, otherwise the deformations are elastic.

The development of the internal state parameter s is described as follows:

$$\dot{\mathbf{s}} = \oplus \mathbf{h}\_{\mathbf{o}} \Big| \mathbf{1} - \frac{\mathbf{s}}{\mathbf{s}^\*} \Big| \stackrel{\mathbf{a}}{\dot{\varepsilon}}\_{\mathbf{F}^{\mathbf{l}}} \tag{7}$$

where a and ho are constants, while s∗ represents the saturated value of the internal parameter. The ⊕ operator is defined to return +1 if s ≤ s∗, otherwise it will return −1. The effect of softening or hardening is included in the model by this operator. The saturation values of <sup>s</sup><sup>∗</sup> depend on the rate of equivalent plastic deformation . ε a pl and can be expressed as follows:

$$\mathbf{s}^\* = \hat{\mathbf{s}} \left\{ \frac{\dot{\boldsymbol{\varepsilon}}^\mathbf{a}}{\mathbf{A}} \mathbf{e}^{(\frac{\mathbf{Q}}{\mathbf{H}})} \right\}^n,\tag{8}$$

where s and n represent constants. ˆ

Expression (7) shows that the development of the parameter s depends on the rate of the equivalent plastic deformation, and at the same time on the current state of the internal state parameter s.

#### **4. Theoretical Description of the Identification Procedure**

The material parameter identification can be described mathematically as:

$$\begin{aligned} \mathsf{f}(\mathsf{X}) &= \text{minimum}, \\ \text{subject to } \mathsf{g}\_{\mathsf{j}}(\mathsf{X}) &> 0, \ \mathsf{j} = \{1, 2, 3, \ldots, \ N\_{\mathsf{C}}\}, \end{aligned} \tag{9}$$

where f(X) is an objective function, X represents a vector of material parameters, gj (X) is j-th constrain function, and NC is the number of constraint functions. For minimalization, the finite element model updating (FEMU) approach was used, which is described in [23]. Schematically, the parameters of the material model were determined by repeated optimization calculations with custom software written in Python programming language. The solution procedure was as follows: FE models of the experiments were created at first. The models were created in commercial software (ANSYS©) as boxes, where the inputs were values of the material parameters and the outputs were values corresponding to loaded or measured data (in this paper these were a force, a displacement, and a current time). The difference between the simulation outputs and the measured data defines the value of an objective function. The simulations were solved in a cycle, whereby the inputs were changed with respect to the minimized value of the objective function. All experiments were deformation-controlled, but forces were used for formulation of the objective function. The difference between experimental data and data from the simulation model for one experiment was solved as an individual objective function fi(X) for the i-th experiment:

$$\mathbf{f}\_{\mathbf{i}}(\mathbf{X}) = \frac{\sum\_{j=1}^{N\_{\mathbf{i}}} \left| \mathbf{F}\_{\mathbf{j}}^{\text{EXP}} - \mathbf{F}\_{\mathbf{j}}^{\text{FEM}}(\mathbf{X}) \right|}{\sum\_{j=1}^{N\_{\mathbf{i}}} \left| \mathbf{F}\_{\mathbf{j}}^{\text{Exp}} \right|}, \text{ i } = \{1, 2, 3, \dots, N\}, \tag{10}$$

where Ni is the number of measurement points for the i-th experiment, FEXP <sup>j</sup> is an experimental force, FFEM <sup>j</sup> (X) is the force obtained from the simulation, and N is the number of experiments.

The value of objective function for all experiments (f(X)) was solved as:

$$\mathbf{f}(\mathbf{X}) = \sqrt{\frac{\sum\_{i=1}^{\mathbf{N}} \mathbf{f}\_i(\mathbf{X})^2}{\mathbf{N}}},\tag{11}$$

A genetic algorithm (GA) was used to identify the material parameters, as in [1]. A GA for the identification of material parameters was used in [24]. The GA was used because it is not dependent on the material model. A detailed description of the GA can be found in [25]. A chromosome is defined by a vector of genes:

$$\mathcal{X} = \{ \mathbf{p}\_{i'} \, \mathbf{i} = \{ 1, 2, 3, \dots, N\_{\mathbb{P}} \} \}, \tag{12}$$

where pi is the i-th gene, which corresponds to the i-th material parameter, and Np is the number of parameters.

An initial population was given by NG solutions, which was generated randomly by using the hill-climbing algorithm [26] from an initial chromosome. For the i-th gene we used:

$$\mathbf{p}\_{\rm i}^{\rm New} = \mathbf{p}\_{\rm i}^{\rm Old} \cdot (1 + \text{Rand} (-\mathbf{k}\_{\rm HC}, \, \mathbf{k}\_{\rm HC})),\\\mathbf{i} = \{1, 2, 3, \dots, N\},\tag{13}$$

where pi Old is the original value of the i-th gene, pi New is the proposed value of the i-th gene, and the Rand function generates random values with uniform distribution from the interval given by ±kHC.

The constraint functions are defined as:

$$\mathbf{g}\_{\mathbf{j}}(\boldsymbol{\lambda}) = \mathbf{p}\_{\mathbf{j}'} \mathbf{j} = \{1, 2, 3, \dots, \text{ N}\_{\mathbb{P}}\},\tag{14}$$

where the number of constraint conditions N*c* is the same as the number of parameters Np.

In each cycle a new individual is created, its objective function value is calculated, and then the individual is added to the population. Regarding the quality of an individual, the solution is determined by the value of the objective function (Equation (11)). The new individual is created in two ways:


The child chromosome was created from three chromosomes of parents, whereby the parents were selected randomly (uniform distribution) from the population. The parents were sorted by the value of the objective function and are indicated by a superscript of 1, 2, or 3; parent 1 had the best (the lowest) objective function and parent 3 had the worst (the greatest) objective function. The child value of the gen (pi Child) is calculated as:

$$\mathbf{p}\_{\text{i}}^{\text{C.hild}} = \mathbf{p}\_{\text{i}}^{\text{Test}} + \boldsymbol{\Delta}\mathbf{p}\_{\text{i}}.\ \mathrm{i} = \{1, 2, 3, \dots, N\},\tag{15}$$

where Δpi is the gene modification value. The modification value is calculated from parent genes as follow:

$$\begin{aligned} \Delta \mathbf{p}\_{\mathbf{i}} &= \mathbf{k} \left( \mathbf{p}\_{\mathbf{i}}^{\mathbf{1}} - \mathbf{p}\_{\mathbf{i}}^{\mathbf{2}} \right), \text{ or } \\ \Delta \mathbf{p}\_{\mathbf{i}} &= \mathbf{k} \left( \mathbf{p}\_{\mathbf{i}}^{\mathbf{1}} - \mathbf{p}\_{\mathbf{i}}^{\mathbf{3}} \right), \text{ i } = \{1, 2, 3, \dots, N\}, \end{aligned} \tag{16}$$

where equations are selected randomly (with 70% and 30% probability); pi 1, pi 2, pi <sup>3</sup> are the gen values of the first parent, second parent, and third parent, respectively; k is a growth coefficient (1 − 2).

The size of the resulting gene values is controlled. If the gene values change too little Δpi < pi Best (<sup>1</sup> ± 0.0001) or too much <sup>Δ</sup>pi > pi Best (<sup>1</sup> ± 0.05), then the new gene values are generated randomly using the hill-climbing algorithm with kHC = 0.001.

Parameters suitable for identification are selected using sensitivity analysis. This sensitivity analysis shows an effect of the parameter on the value of the objective function with respect to the other parameters. The sensitivity analysis is performed from Np + 1 calculations as follows:

$$\mathbf{f} = \mathbf{f}(\mathbf{X}),$$

$$\mathbf{f}\_{+\Delta \mathbf{i}} = \mathbf{f}\left(\left\{\mathbf{p}\_{1'}\mathbf{p}\_{2'}, \dots, \mathbf{p}\_{i}(1+\mathbf{k}\_{\text{Sen}}), \dots, \mathbf{p}\_{\text{Np}}\right\}\right),\tag{17}$$

where kSen is a sensitivity coefficient (0.05 − 0.001). The sensitivity value Si is calculated according to the equation:

$$\text{Si} = \frac{100}{\text{k}\_{\text{Sen}}} \left| 1 - \frac{\text{f}\_{+\text{Ai}}}{\text{f}} \right| \left[ \text{\textdegree\text{\textdegree\textdegree\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text}} \left[ \text{\textquotestext{\textquotestext{\textquotestext{\textquotesbltext}}}{} \text{\textquotesingle}} \right] \left[ \text{\textquotesingle\text\textquotesingle} \right] \text{\textquotesingle} = \{1, 2, 3, \dots, \text{N} \text{\textquotesingle} \text{\textquotesingle}\}. \tag{18}$$

Parameters whose sensitivity values are greater than the critical value SKrit are selected for identification:

$$\mathbf{S}\_{\mathbf{i}} \rhd \mathbf{S}\_{\mathbf{K} \text{rit}} \rightarrow \text{identity} \\ \mathbf{p}\_{\mathbf{i}\prime} \mathbf{i} = \{1, 2, 3, \dots, N \mathbf{p}\} . \tag{19}$$

In later phases of identification, it is advisable to check if the given parameter pi is not in the local minimum of the objective function. Here, the number of required calculations is 2Np + 1:

$$\mathbf{f}\_{+\Lambda\dot{\mathbf{i}}} > \mathbf{f} < \mathbf{f}\_{-\Lambda\dot{\mathbf{i}}\prime} \ \mathbf{i} = \{1, 2, 3, \dots, N\mathbf{p}\}.\tag{20}$$

#### **5. Description of the Identification Process**

Ten experiments were used for identification (see Figure 3a–f), however the experiments at 44 ◦C were excluded from the identification process and were used only for validation. The assignment of the experiments to the values of the individual objective functions is shown in Table 1.

**Table 1.** Assignment of the individual objective functions to individual experiments.


In Table 1, the temperatures are denoted as T1 = 23 ◦C, T2 = 44 ◦C, T3 = 60 ◦C, and T4 = 80 ◦C; the deformation rates are . x1 <sup>=</sup> 0.017 mm s<sup>−</sup>1, . x2 = 0.167 mm s<sup>−</sup>1, and . x3 = 1.667 mm s<sup>−</sup>1; ID denotes the indentation tests and GT denotes the graded tensile tests (see Chapter 2). Data from the tests at temperature T2 were used to validate the parameters.

The initial values of parameters are shown in Table 2.

**Table 2.** Initial parameters of the Anand material model.


The value of the modulus of elasticity E is temperature-dependent, which is shown in Figure 3a. For this reason, three values of the modulus of elasticity E20, E50, and E80 (for 20, 50, and 80 ◦C) were added to the material model. The temperature chamber made it impossible to use optical methods (e.g., the digital image correlation method (DIC)) [27] to measure the Poisson number, so this parameter was excluded from identification.

Chromosome (12) consists of these material parameters, whereby the individual genes correspond to the individual material parameters:

$$\mathcal{X} = \left\{ \mathbf{E}\_{20\prime} \mathbf{E}\_{50\prime} \mathbf{E}\_{80\prime} \mathbf{s}\_{0\prime} \frac{\mathbf{Q}}{\mathbf{R}}, \mathbf{A}\_{\prime} \mathbf{x}\_{\mathrm{i}\prime} \mathbf{m}\_{\prime} \mathbf{h}\_{0\prime} \mathbf{\hat{S}}\_{\prime} \mathbf{n}\_{\prime} \mathbf{a} \right\}. \tag{21}$$

The hill-climbing algorithm was used with the parameter kHC = 0.001, the genetic algorithm was used with the parameter NG = 10, and the growth coefficient k = 1.5 was used in both cases. The sensitivity coefficient was kSen = 0.01. The sensitivity value shown in (18) and the criterion shown in (20) were recalculated. The way in which the model parameters were determined is described by the following three tables. Table 3 contains the parameter values during the solution.

**Table 3.** Parameter values with three moduli of elasticity during identification.


The values for the objective function during the solution are shown in Table 4 and the sensitivity values are shown in Table 5. The initial parameters for identification (values of the objective function) are given in line 0 of Table 3. The sensitivity was calculated from "line 0" values, and the resulting values are shown on line 1 of Table 5. The genetic algorithm was run 200 times in every identification cycle. This process was repeated in lines 2 to 5. Line 5 in Tables 4 and 5 shows the resulting parameter and error values. The fulfilment of criterion (20) is indicated in Table 5 as U. The resulting values (line 5, Table 4) were used to calculate the sensitivity (line 6\* in Table 5). Values in line 6\* were no longer used for further calculation. This shows that the minimum given by the coefficient kSen had not been reached yet. Table 5 shows in the last step (line 5) a decrease of the objective function value by 0.1% after 200 calculation cycles; the identification process was, therefore, completed.

**Table 4.** The objective function values during identification.



**Table 5.** Sensitivity values during identification. Parameters that fulfil the criterion are shown in bold.

#### **6. Results**

The obtained final parameter values of the Anand material model are summarized in Table 6 and the FEM analyses were performed with these material parameters. The value of the Poisson ratio was *μ* = 0.33. Figures 4 and 5 show a comparison of experiments with results from the FEM analysis.

**Table 6.** Final parameter values.


**Figure 4.** Comparison of experiments with FEM solutions: (**a**) indentation test at 23 ◦C; (**b**) tensile test at 23 ◦C and rate of deformation of 1.667 mm s<sup>−</sup>1.

**Figure 5.** *Cont*.

**Figure 5.** Comparison of experiments with FEM solutions (continuation): (**a**) tensile test at 60 ◦C and rate of deformation of 0.017 mm s−1; (**b**) tensile test 60 ◦C and rate of deformation of 0.167 mm s−1; (**c**) indentation test at 60 ◦C; (**d**) tensile test at 60 ◦C and rate of deformation of 1.667 mm s−1; (**e**) graded tensile test at 60 ◦C; (**f**) tensile test at 80 ◦C and rate of deformation of 0.017 mm s<sup>−</sup>1; (**g**) indentation test at 80 ◦C; (**h**) graded tensile test at 80 ◦C.

Very good agreement can be seen between the experiments and the FEM solution. Different behavior just before the failure of the specimen is shown in Figure 5a,b,d–f,h. This behavior is less significant under lower temperatures and higher deformation rates. The difference can be seen in Figure 5a,f especially (temperatures of 60 and 80 ◦C and deformation rate of 0.017 m s<sup>−</sup>1). This effect was not captured in the simulations.

#### **7. Validation**

Validation testing was carried out on a series of experiments at 44 ◦C. Experiments at this temperature were not used for material parameter identification. The results are shown in Figure 6. In Figure 6, a very good agreement can be seen between FEM solutions and experiments. The greater error occurs just before the failure of the specimen, as also mentioned in the previous chapter. The conclusion is that the material models with those material parameters do not accurately describe the mechanical behavior shortly before the sample ruptures. In order to describe this area of difference better, more experiments will be necessary, especially focused on this small area. If no experiments of this kind are carried out, the influence of this small area on the value of the objective function is negligible.

**Figure 6.** Comparison of FEM solutions with validation experiments at 44 ◦C: (**a**) simple tensile test at rate of deformation of 0.017 mm s−1; (**b**) simple tensile test at rate of deformation of 0.167 mm s−1; (**c**) graded tensile test; (**d**) indentation experiment.

#### **8. Discussion**

The GA modification was used to identify the parameter values. In order to reduce the number of calculation cycles and individuals in the GA population, the number of parameters was reduced by means of sensitivity control. The sensitivity calculation was performed at the beginning of the solution, and then again when the GA convergence rate decreases significantly.

Further sensitivity tests showed that if some parameters are changed there is still a decrease, while others are at their minima in the range tested (see table of resulting sensitivity U). In subsequent cycles, the value of the objective function gradually decreased in the range of 0.01% to 100 GA steps. This value decreased further in each cycle. The resulting effects of these calculations on the value of the objective function and the values of parameters were negligible, and therefore are not described in this contribution.

The Anand constitutive model used in this article contained 9 parameters and 3 moduli of elasticity, which were identified for temperatures of 20, 50, and 80 ◦C. Figure 3a shows the influence of temperatures on the value of E. This dependence appeared at all speeds tested, so it was included in the identification, although its effect was not significant. Input values of the material parameters were taken from [1] for the ABS-M30 material and determined at a temperature of 23 ◦C and displacement rate of 1 mm s−<sup>1</sup> from tensile tests and indentation tests, while the Poisson's ratio was determined from tensile tests with DIC. These experiments were not used for identification in this contribution because they showed large differences (originally E = 1780 MPa, new value E = 1100 MPa) in the behavior compared to the results presented above. The authors believe that this was due to a difference in the supply of the material and different 3D printing settings.

In [1], indentation and DIC measurement data were also used to identify and validate the results. When using a temperature chamber, it was not possible to perform DIC measurements—more extensive modifications of the experimental equipment would be necessary. In the case of DIC, this mainly involves determining the value of the Poisson ratio, which was, therefore, taken from [1]. Indentation tests describe the behavior of multiaxial stresses in the area above the yield point in the compression area. Tensile tests carried out at lower temperatures and higher speeds (Figure 3b,c) show almost linear behavior.

The use of different types of experiments (e.g., tensile tests, indentation tests) for identification or at least for validation is the key for the further use of the identified parameters. Indentation, tensile, and graded tensile tests at 44 ◦C were used for validation (see [28]). Components produced using 3D printing are not only loaded by tension, so the authors assume that after modifying the experimental equipment, a different set of experiments will be performed. On the other hand, conventional components are loaded in the area corresponding to the first half of the loading curve, where tensile tests in particular show very good agreement with the FEM solution.

#### **9. Conclusions**

The obtained results show the possibility of finding the material model parameters using indentation or tensile tests. The results of the FEM calculations are in good agreement with the data obtained experimentally. Attention will be paid to the modification of the test equipment in order to enable DIC measurements. The Anand material model can be used for specimens loaded at different temperatures. It was found that the Anand material model describes the behavior of the ABS-M30 material very well with respect to different temperatures and loading speeds. On the other hand, the Anand model is not able to capture the anisotropic behavior of the investigated material sufficiently. For tests on specimens printed in the most unfavorable position (related to the load capacity of the specimen), the material model parameters identified as described above will give conservative results and the model can be used to design components. Studies on a more advanced material model that includes the effects of material anisotropy and a damage model are currently in progress. The proposed identification methodology is not dependent on the material model.

**Author Contributions:** Conceptualization, Z.P.; material model description, M.F.; planning and execution of experiments, F.F.; identification procedure and FEM, J.R.; resources, M.F., J.R., and F.F.; writing—original draft preparation, Z.P., M.F., and J.R.; writing—review and editing, Z.P., M.F., J.R., and F.F.; visualization, Z.P. and F.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data can be provided upon request from the correspondent author.

**Acknowledgments:** This article was elaborated under support of the Research Center of Advanced Mechatronic Systems, reg. no. CZ.02.1.01/0.0/0.0/16\_019/0000867, within the framework of the Operational Program for Research, Development, and Education.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

