*Article* **Hardening Parameter Homogenization for J2 Flow with Isotropic Hardening of Steel Fiber-Reinforced Concrete Composites**

**Petr V. Sivtsev 1,\* and Piotr Smarzewski <sup>2</sup>**


**Abstract:** Numerical modeling of the stress–strain state of composite materials such as fiber-reinforced concrete is a considerable computational challenge. Even if a computational grid with the resolution of all inclusions is built, it will take a great amount of time for the most powerful clusters to calculate the deformations of one concrete block with ideal parallelization. To solve this problem, the method of numerical homogenization is actively used. However, when plastic deformations are taken into account, the numerical homogenization becomes much more complicated due to nonlinearity. In this work, the description of the anisotropic nature of the hardening of the composite material and the numerical homogenization for the J2 flow with isotropic hardening is proposed. Here, the deformation of a composite material with a periodic arrangement of inclusions in the form of fibers is considered as a model problem. In this case, the assumption is made that inclusions have pure elastic properties. Numerical homogenization of the elasticity and plasticity parameters is performed on the representative element. The novelty of the work is related to the attempt at hardening parameter homogenization. The calculated effective parameters are used to solve the problem on a coarse mesh. The accuracy of using the computational algorithm is checked on model problems in comparison with the hardening parameters of the base material. The finite element implementation is built using the FEniCS computing platform and the fenics-solid library.

**Keywords:** elastoplastic; hardening; homogenization; J2 flow; composites; fiber reinforcement; finite element

#### **1. Introduction**

Composite materials are used in various industries. Over time, their structures change, and their production is optimized in terms of operational properties. Composite materials have the advantage of a high strength-to-weight ratio, improved thermal conductivity, or isolation. Most composites are materials with elastoplastic properties, for example, reinforced concrete [1–3]. The behavior of traditional materials has been studied well enough, and there are mathematical models verified over time that describe the deformation properties of a particular natural material, whereas the description of composite materials is not so simple. In particular, the description of the deformations of a composite material outside the theory of linear elasticity can be handled by a second-order radial return algorithm [4].

In our study, we consider metal inclusions. For the mathematical description of the deformation of concrete, we use the theory of plastic flow with isotropic hardening and a description of the plastic flow region using the Mises yield criterion. Numerical modeling of a composite material with a mesh resolution of all inclusions in the form of reinforcement and distributed fibers is a considerable computational problem [5,6]. Modern computer

**Citation:** Sivtsev, P.V.; Smarzewski, P. Hardening Parameter Homogenization for J2 Flow with Isotropic Hardening of Steel Fiber-Reinforced Concrete Composites. *Crystals* **2021**, *11*, 776. https://doi.org/10.3390/cryst11070776

Academic Editor: Tomasz Sadowski

Received: 8 June 2021 Accepted: 30 June 2021 Published: 2 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

equipment has impressive power due to a large number of cores with a high clock speed. However, the complexity of the problem increases nonlinearly with an increase in the structure size under consideration.

Previous work [7] proposed the use of the effective tensor of elasticity, obtained from numerical homogenization, the parameters of plasticity, such as the yield stress, and the hardening parameter as in concrete. At small values of plastic deformations, this approximation demonstrates good accuracy and can be used to describe plastic deformations until the body begins to collapse. In other words, our model considers concrete stress–strain dependence as piecewise linear, as shown in Figure 1. This model is supposed to be used for stresses below material strength [8].

**Figure 1.** Plastic flow with isotropic hardening representation (blue).

In this work, a description of the anisotropic nature of the hardening of the composite material and numerically homogenizing the hardening parameter for the J2 flow with isotropic hardening is proposed. Deformation of a composite material with a periodic arrangement of inclusions in the form of fibers is considered as a model problem.

On the representative element, the elastic parameters and the hardening parameter are numerically homogenized. The calculated effective parameters are used to solve the problem on a coarse mesh. The accuracy of the computational algorithm is checked on model problems in comparison with the use of the hardening parameter of the base material, i.e., concrete. The finite element implementation is built using the FEniCS [9] computing platform and the fenics-solid library [10].

#### **2. Materials and Methods**

#### *2.1. Problem Statement*

A two-dimensional mathematical model describing the stress–strain state of a concrete composite with steel fibers is considered, described by the equilibrium equation

$$\text{div}\sigma = 0,\ \text{x} \in \Omega = \Omega\_1 \cup \Omega\_2. \tag{1}$$

where Ω<sup>1</sup> is a concrete subdomain, Ω<sup>2</sup> is a subdomain of fiber inclusions, *σ* is a stress tensor. For the plasticity model [11,12], the elastic tensor consists of elastic and plastic parts:

$$\mathfrak{e} = \begin{pmatrix} \frac{\partial u\_1}{\partial x\_1} \\ \frac{\partial u\_2}{\partial x\_2} \\ \frac{\partial u\_2}{\partial x\_1} + \frac{\partial u\_1}{\partial x\_2} \end{pmatrix} = \mathfrak{e}^c + \mathfrak{e}^p. \tag{2}$$

Next, the connection between the stress tensor and the elastic strain tensor using the generalized Hooke's law is added. For convenience, Voight notation is used. In a two-dimensional case, Hooke's law can be presented as:

$$
\sigma = \begin{pmatrix} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{12} \end{pmatrix} = \mathbb{C}\mathfrak{E}^{\mathfrak{E}}, \quad \mathbb{C} = \begin{pmatrix} \mathbb{C}\_{1111} \, \mathbb{C}\_{1122} \, \mathbb{C}\_{1112} \\ \mathbb{C}\_{2211} \, \mathbb{C}\_{2222} \, \mathbb{C}\_{2212} \\ \mathbb{C}\_{1211} \, \mathbb{C}\_{1222} \, \mathbb{C}\_{1212} \end{pmatrix}, \tag{3}
$$

where *C* is an elasticity tensor, which is as follows for isotropic materials:

$$\mathbf{C} = \begin{pmatrix} \lambda + 2\mu & \lambda & \mathbf{0} \\ \lambda & \lambda + 2\mu & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mu \end{pmatrix},\tag{4}$$

where *λ* and *μ* are the Lamé coefficients determined through Young's modulus *E* and Poisson's ratio *ν*. Model values of these parameters for concrete and steel are presented in Table 1 [13,14].

**Table 1.** Model elastic parameters of concrete matrix, basalt, and steel fibers.


Yield criterion of plastic flow with isotropic hardening is following:

$$f(\sigma, \,\varepsilon^p, \,\kappa) = \phi(\sigma) - q\_{\text{iso}}(\kappa) - \sigma\_{\text{y}} \le 0,\tag{5}$$

where *φ*(*σ*) is a scalar effective stress measure, *qiso*(*κ*) is a scalar stress-like internal variable used to model isotropic hardening, *κ* is an internal variable, and *σ<sup>y</sup>* is initial yield stress. In this work, a von Mises model (also known as J2 flow) with linear isotropic hardening is used. For this model, *φ* and *qiso* are written as:

$$\phi(\sigma) = \sqrt{\frac{3}{2} s\_{ij} s\_{ij\prime}}\tag{6}$$

$$q\_{iso}(\kappa) = H\kappa,\tag{7}$$

where *sij* = *σij* − *σkkδij*/3 is the deviatoric stress and the constant *H* is a hardening parameter. For associative von Mises plasticity with isotropic hardening, there are the following rates:

$$
\dot{\varepsilon}^p = \dot{\lambda} \frac{\partial f}{\partial \sigma'} \tag{8}
$$

$$
\kappa = \sqrt{\frac{2}{3} \dot{\varepsilon}\_{ij}^p \dot{\varepsilon}\_{ij}^p} \tag{9}
$$

Concrete hardening parameters *H* and yield stress σ<sup>y</sup> are taken as equal to 21.5 GPa and 30.0 MPa, respectively. These values are related to compression and used for all kinds of deformation in the modeling purposes of this study.

#### *2.2. Hardening Parameter Homogenization Algorithm*

#### 2.2.1. Effective Hardening Parameter

The main assumption of this work is that the effective hardening parameter depends on the tensor of plastic deformations and takes a certain numerical value for each of its components. In this case, for the classical model of isotropic hardening, the description of compression and tension is described in the same way. The simplest form of the effective hardening factor satisfying these criteria will be as follows:

$$\begin{array}{c|c|c|c|c|c} \mathbf{h}\_{eff} = \frac{|\boldsymbol{\varepsilon}\_{11}^{p}|h\_{1} + |\boldsymbol{\varepsilon}\_{22}^{p}|h\_{2} + |\boldsymbol{\varepsilon}\_{12}^{p}|h\_{3}}{|\boldsymbol{\varepsilon}\_{11}^{p}| + |\boldsymbol{\varepsilon}\_{22}^{p}| + |\boldsymbol{\varepsilon}\_{12}^{p}|} \text{ for } \left|\boldsymbol{\varepsilon}\_{11}^{p}\right| + \left|\boldsymbol{\varepsilon}\_{22}^{p}\right| + \left|\boldsymbol{\varepsilon}\_{12}^{p}\right| \neq 0, \\\ h\_{eff} = h\_{c} \text{ for } \left|\boldsymbol{\varepsilon}\_{11}^{p}\right| + \left|\boldsymbol{\varepsilon}\_{22}^{p}\right| + \left|\boldsymbol{\varepsilon}\_{12}^{p}\right| = 0, \end{array} \tag{10}$$

the following is true for it:

$$\begin{array}{l} \text{h}\_{\varepsilon ff} = \text{h}\_{1} \text{ for } \varepsilon^{p}\_{22} = \varepsilon^{p}\_{12} = 0 \text{ and } \varepsilon^{p}\_{11} \neq 0, \\\text{h}\_{\varepsilon ff} = \text{h}\_{2} \text{ for } \varepsilon^{p}\_{11} = \varepsilon^{p}\_{12} = 0 \text{ and } \varepsilon^{p}\_{22} \neq 0, \\\text{h}\_{\varepsilon ff} = \text{h}\_{3} \text{ for } \varepsilon^{p}\_{11} = \varepsilon^{p}\_{22} = 0 \text{ and } \varepsilon^{p}\_{12} \neq 0, \end{array} \tag{11}$$

Let us take the effective value of the yield point equal to the concrete yield point. Next, consider the algorithm with which one can calculate the values of *h*1, *h*2, *h*3.

#### 2.2.2. Elastic Parameter Homogenization

First, the representative elements are selected. For the unidirectional location of fiber, a unit cell is used. For the bidirectional location of fibers, a 2 × 2 cell is used. In the case with a random distribution of fibers, the accuracy of the elasticity problem highly depends on representative volume size. Therefore, an optimal representative volume with 16 fibers is chosen. The effective elastic tensor of the composite material *Ceff* is calculated using representative volume, as in previous work [15].

#### 2.2.3. Hardening Parameter Homogenization

Second, three problems of elastic plasticity are solved using very large values of the boundary condition, for example, when the deformations on whole boundary are equal:

$$\begin{array}{l} u\_D = \left(10^6 \cdot x\_{1\prime} \cdot 0\right), \\ u\_D = \left(0, \ 10^6 \cdot x\_2\right), \\ u\_D = \left(10^6 \cdot x\_2 / 2, \ 10^6 \cdot x\_1 / 2\right). \end{array} \tag{12}$$

Additionally, for each case, we obtain the averaged stress tensors *σ*1, *σ*2, *σ*3, which must correspond to solutions with effective parameters. Given the large value of deformation and, accordingly, the large value of stresses, it can be assumed that purely elastic deformations are negligible, and the ratio of stress and strain is described by the elasticplasticity tensor. Thus, by solving the same problems on a homogeneous medium using the effective parameters, we must obtain the same stress values. However, since complete correspondence to the stress tensor by varying the value of the hardening parameter alone is hardly obtained, the agreement that it is necessary to obtain correspondence with the most essential components of the tensors is taken. Namely, *σ*<sup>1</sup> <sup>11</sup>, *<sup>σ</sup>*<sup>2</sup> <sup>22</sup>, *<sup>σ</sup>*<sup>3</sup> 12.

Further, by solving the problem on a coarse mesh with a uniform distribution of the effective parameters *Ceff* , *heff* with boundary conditions (12), correspondence with *σ*1 <sup>11</sup>, *<sup>σ</sup>*<sup>2</sup> <sup>22</sup>, *<sup>σ</sup>*<sup>3</sup> <sup>12</sup> must be obtained. When solving the problem with boundary conditions (12), smallness of the elastic deformations occurs and, accordingly, the following becomes true:

$$
\sigma \approx \mathbb{C}^{cp} \varepsilon \tag{13}
$$

At the same time, from the theory of elastic flow with isotropic hardening, the elastoplastic tangent is equal to

$$\mathbf{C}^{\sigma p} = \boldsymbol{\Xi} - \frac{\boldsymbol{\Xi} : \partial\_{\sigma} \boldsymbol{g} \otimes \partial\_{\sigma} f : \boldsymbol{\Xi}}{\partial\_{\sigma} f : \boldsymbol{\Xi} : \partial\_{\sigma} \boldsymbol{g} + h},\tag{14}$$

where Ξ = *C*/ *I* + Δλ*C* : *∂*<sup>2</sup> *σσg* .

As long as the values of the strain tensor for each point are the same, a hardening parameter that will correspond to the selected direction of deformation can be calculated at any point. To calculate the parameters, the following relations can be used:

$$h\_1 = \frac{R^1[1, 1]}{\Xi^1[1, 1] - \frac{\sigma^1\_{11}}{\varepsilon^1\_{11}}} - a^1, \ h\_2 = \frac{R^2[2, 2]}{\Xi^2[2, 2] - \frac{\sigma^2\_{22}}{\varepsilon^2\_{22}}} - a^2, \ h\_3 = \frac{R^3[3, 3]}{\Xi^3[3, 3] - \frac{\sigma^3\_{12}}{\varepsilon^3\_{12}}} - a^3,\tag{15}$$

where *R* = Ξ : *∂σg* ⊗ *∂σ f* : Ξ, and *α* = *∂σ f* : Ξ : *∂σg* (top index corresponds to each problem). Thus, by solving three problems with different boundary conditions (12), the sought-for hardening parameter values can be obtained.

#### *2.3. Research Object*

A two-dimensional concrete structure with inclusions in the form of fibers is considered for different cases of fibers:


The aspect ratio of fibers is 1 to 10 and the concentration is equal to 2.5%. The locations of the fibers for each case are shown in Figure 2. The number of fibers is 256.

**Figure 2.** Research object geometries: (**a**) unidirectional location of fibers; (**b**) bidirectional location of fibers; (**c**) random location of fibers.

#### **3. Results**

For the object of study, the values of the elastic tensor coefficients and the components of the anisotropic hardening coefficient in the case of unidirectional, bidirectional, and random steel fiber distribution were calculated.

To check the effective parameters, the following tests were performed:


The purpose of the tests was to identify the scope of the proposed model.

In the first test, the left border is completely fixed, i.e., a Dirichlet boundary condition equal to zero vector is applied, and normal stress is applied to the right border, which increases from 0 to 200 MPa over 100 time steps. The second test is similar to the first, but anchoring is at the bottom, and tensile stress is on the top. The third test is related to tangential stress at the right side with anchoring on the left side.

In contrast to other tests, the fourth test features compression along two axes. A quarter of the body with 1024 fibers is considered. For this problem, the random distribution is modeled as symmetric for vertical and horizontal axes. Thus, the simulated geometry is the same. Normal stress is set above and on the right, and it also increases from 0 to 200 MPa over 100 time steps. The corresponding symmetry boundary conditions are set on the left and below.

**Figure 3.** Schematic description of the tests: (**a**) Schematic description of test 1. Full anchoring of the left side and normal compression from the right. (**b**) Schematic description of test 2. Full anchoring from the bottom edge and normal compression from the top.

**Figure 4.** Schematic description of the tests: (**a**) Schematic description of test 3. Full anchoring of the left side and tangential compression from the right. (**b**) Schematic description of test 4. Symmetry conditions on the left and bottom faces and normal tension on the right and top faces.

In addition to material variation, the solution of the elastic plasticity problem with isotropic hardening is performed on an inhomogeneous region on a computational grid with full geometry (Figure 5a–c). This solution is considered as an exact solution. The errors of the compared models are calculated referring to it. A coarse grid (Figure 5d) requires a solution made using the effective elastic tensor. The models using the effective hardening factor and the concrete hardening parameter are compared. When using return mapping, the value of the plastic strain tensor for calculating the effective coefficient is taken from the previous iteration.

**Figure 5.** Computational grids: (**a**) Computational grid for case with unidirectional fiber distribution. (**b**) Computational grid for case with bidirectional fiber distribution. (**c**) Computational grid for case with random fiber distribution. (**d**) A coarse computational grid on which the numerical solution was performed using the effective parameters.

The computational mesh of the complete geometry for unidirectional, bidirectional, and random fiber distribution and the coarse computational mesh, on which the solution is performed using the effective parameters, have 33,963 vertices, 78,361 cells, 44,814 vertices, 100,063 cells, 47,822 vertices, 106,910 cells, and 2319 vertices, 4640 cells, respectively.

#### *3.1. Calculation of Effective Parameters*

To calculate the effective parameters, for the unidirectional location of fibers, a unit cell is used. For the bidirectional location of fibers, we used a 2 × 2 cell. In the case with a random distribution of fibers, the accuracy of the elasticity problem highly depends on the representative volume size. It is thought that a 16-fiber case would be optimal. The computational meshes for calculating the elastic parameters are shown in Figure 6a–f and the computational mesh for calculating the hardening parameters is shown in Figure 6c.

The computational mesh for calculating the elastic parameters for unidirectional, bidirectional, and random distribution and the coarse computational mesh for calculating the hardening parameters have 528 vertices, 1126 cells, 2020 vertices 4319 cells, 2627 vertices, 5960 cells and 177 vertices, 312 cells, respectively.

After applying the algorithm, the effective parameters presented in Table 2 were obtained.

**Figure 6.** Computational grids used to calculate the effective parameters: (**a**) Computational mesh of a representative element for unidirectional case. (**b**) Computational mesh of a representative element for bidirectional case. (**c**) Coarse computational mesh used to calculate the effective hardening parameters. (**d**) Computational mesh of a representative element with random distribution 1. (**e**) Computational mesh of a representative element with random distribution 2. (**f**) Computational mesh of a representative element with random distribution 3.


**Table 2.** Calculated effective elasticity and plasticity parameters for all cases: unidirectional, bidirectional, and random distribution of steel fibers.

#### *3.2. Compression along the Horizontal Axis*

Based on the results of test 1, the dependences of the maximum value of the displacement modulus in all cases were obtained. They are presented in Figures 7–9.

**Figure 7.** A graph of the dependence of the maximum value of the deformation modulus in Table 1 for unidirectional fiber location case.

**Figure 8.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 1 for bidirectional fiber location case.

The distributions of the deformation error modulus |*uerr*|, which can be expressed by fine mesh solution *ufine* taken as exact and homogenized solution *uhom* as <sup>|</sup>*uerr*<sup>|</sup> <sup>=</sup> *uX fine* − *<sup>u</sup><sup>X</sup> hom*<sup>2</sup> + *uY fine* − *<sup>u</sup><sup>Y</sup> hom*<sup>2</sup> , for test 1 for unidirectional, bidirectional, and random distribution of fibers, are shown in Figures 10–14, respectively.

**Figure 10.** Distribution of the deformation error modulus for test 1 for unidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 11.** Distribution of the deformation error modulus for test 1 for bidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 12.** Distribution of the deformation error modulus for test 1 for random distribution case 1: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 13.** Distribution of the deformation error modulus for test 1 for random distribution case 2: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 14.** Distribution of the deformation error modulus for test 1 for random distribution case 3: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

In all cases, it can be seen that when using the effective hardening parameter, the solution has a significantly lower value of the deformation error.

#### *3.3. Compression along the Vertical Axis*

Based on the results of test 2, the dependences of the maximum value of the displacement modulus in all cases were obtained. They are presented in Figures 15–17.

**Figure 15.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 2 for unidirectional fiber location case.

**Figure 16.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 2 for bidirectional fiber location case.

**Figure 17.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 2 for random fiber distribution case.

The distributions of the deformation modulus error for test 2 for unidirectional, bidirectional, and random distribution of fibers are shown in Figures 18–22, respectively.

**Figure 18.** Distribution of the deformation error modulus for test 2 for unidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 19.** Distribution of the deformation error modulus for test 2 for bidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 20.** Distribution of the deformation error modulus for test 2 for random distribution case 1: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 21.** Distribution of the deformation error modulus for test 2 for random distribution case 2: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 22.** Distribution of the deformation error modulus for test 2 for random distribution case 3: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

In both unidirectional and bidirectional cases, as in test 1, it can be seen that the method proposed in the article gives a solution with a significantly smaller error value. However, for the case with random distribution error values when using effective and base hardening, parameters are similar. This could be due to a bad representative volume.

#### *3.4. Tangential Stress*

Based on the results of test 3, the dependences of the maximum value of the displacement modulus in all cases were obtained. They are presented in Figures 23–25.

**Figure 23.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 3 for unidirectional fiber location case.

**Figure 24.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 3 for bidirectional fiber location case.

The distributions of the deformation modulus error for test 3 for unidirectional, bidirectional, and random distribution of fibers are shown in Figures 26–30, respectively.

**Figure 26.** Distribution of the deformation error modulus for test 3 for unidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 27.** Distribution of the deformation error modulus for test 3 for bidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 28.** Distribution of the deformation error modulus for test 3 for random distribution case 1: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 29.** Distribution of the deformation error modulus for test 3 for random distribution case 2: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 30.** Distribution of the deformation error modulus for test 3 for random distribution case 3: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

In all cases, as in previous tests, we see that the method proposed in the article gives a solution with better accuracy. In particular, this can be seen for random distribution. However, accuracy could be worse for other representative volumes.

#### *3.5. All-Round Compression*

Based on the results of test 4, the dependences of the maximum value of the displacement modulus in all cases were obtained. They are presented in Figures 31–33.

**Figure 31.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 4 for unidirectional fiber location case.

**Figure 32.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 4 for bidirectional fiber location case.

**Figure 33.** A graph of the dependence of the maximum value of the deformation modulus on the traction value for solving the problem of plasticity, taking into account inclusions, using effective parameters, effective parameters of elasticity, and plastic flow parameters of concrete for test 4 for random fiber distribution case.

The distributions of the deformation modulus error for test 4 for unidirectional,

**Figure 34.** Distribution of the deformation error modulus for test 4 for unidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 35.** Distribution of the deformation error modulus for test 4 for bidirectional case: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 36.** Distribution of the deformation error modulus for test 4 for random distribution case 1: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 37.** Distribution of the deformation error modulus for test 4 for random distribution case 2: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

**Figure 38.** Distribution of the deformation error modulus for test 4 for random distribution case 3: (**a**) When using effective hardening parameters; (**b**) when using the concrete hardening parameter.

In contrast to all other tests, the results of the fourth test show the weak side of the proposed model, namely, when comparable values of the strain tensor components occur. In this case, it is impossible to take into account the whole anisotropic nature of the elastic–plastic flow of the composite material.

#### **4. Discussion**

J2 flow with isotropic hardening can be used to model plastic behavior of a great variety of materials, for example, concrete. Fiber-reinforced concrete is currently of great interest. Additionally, to model building blocks made of fiber-reinforced concrete, proper homogenization or multiscale methods are required. The main focus of such methods is to homogenize the elastoplastic tangent tensor. The algorithm presented here of hardening coefficient homogenization is simple to understand in comparison to other plasticity problem homogenization algorithms.

In a previous study, the algorithm of the homogenization of fiber-reinforced concrete with elasticity tensor homogenization using a hardening parameter and the yield stress of concrete without reinforcement was investigated. A satisfactory result with good accuracy was obtained. In this study, we tried to add hardening parameter homogenization in order to improve accuracy. A comparison with the previous algorithm was performed using three different fiber locations and four tests.

#### **5. Conclusions**

According to the results obtained, the proposed model of the hardening coefficient homogenization demonstrates satisfactory results when one of the components of the strain tensor prevails. Under this condition, the results obtained have an error less than when using the concrete hardening coefficient, as shown by tests 1, 2, and 3.

However, when the strain tensor components have comparable values, error values become similar to the values of error of the simple model from previous work. This is due to the anisotropy of the plastic flow that cannot be fully taken into account in a simple numerical change in the hardening coefficient. In this case, we note that both approximations work quite accurately at small plastic deformations.

In further work, it will be necessary to obtain a method of numerical homogenization, which will have more degrees of freedom and describe the anisotropy of plastic deformations better. Different values of plastic parameters of concrete, for compression and tension, should also be addressed in future work.

**Author Contributions:** Conceptualization, P.V.S. and P.S.; Data curation, P.V.S.; Formal analysis, P.V.S. and P.S.; Funding acquisition, P.S.; Investigation, P.V.S. and P.S.; Methodology, P.V.S. and P.S.; Software, P.V.S.; Supervision, P.S.; Validation, P.V.S. and P.S.; Visualization, P.V.S.; Writing—original draft, P.V.S.; Writing—review and editing, P.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by a grant of the President of the Russian Federation, grant number MK-2249.2021.1.1 and "Subvention for Science" (MEiN), project no. FN-4/2021.

**Acknowledgments:** The authors would like to thank every person/department who helped thorough out the research work. The careful review and constructive suggestions by the anonymous reviewers were gratefully acknowledged.

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

#### **References**

