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

: Numerical modeling of the stress–strain state of composite materials such as ﬁber-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 ﬂow with isotropic hardening is proposed. Here, the deformation of a composite material with a periodic arrangement of inclusions in the form of ﬁbers 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 ﬁnite element implementation is built using the FEniCS computing platform and the fenics-solid library.


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][2][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 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].

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 where Ω is a concrete subdomain, Ω 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: 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].

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 where Ω 1 is a concrete subdomain, Ω 2 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: Crystals 2021, 11, 776 3 of 26 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: where C is an elasticity tensor, which is as follows for isotropic materials: 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]. Yield criterion of plastic flow with isotropic hardening is following: where φ(σ) is a scalar effective stress measure, q iso (κ) is a scalar stress-like internal variable used to model isotropic hardening, κ is an internal variable, and σ y is initial yield stress. In this work, a von Mises model (also known as J 2 flow) with linear isotropic hardening is used. For this model, φ and q iso are written as: where s ij = σ 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: Concrete hardening parameters H and yield stress σ y 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.

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 Crystals 2021, 11, 776 4 of 26 compression and tension is described in the same way. The simplest form of the effective hardening factor satisfying these criteria will be as follows: the following is true for it: 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 .

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 C e f f is calculated using representative volume, as in previous work [15].
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, σ 1 11 , σ 2 22 , σ 3 12 . Further, by solving the problem on a coarse mesh with a uniform distribution of the effective parameters C e f f , h e f f with boundary conditions (12), correspondence with σ 1 11 , σ 2 22 , σ 3 12 must be obtained. When solving the problem with boundary conditions (12), smallness of the elastic deformations occurs and, accordingly, the following becomes true: At the same time, from the theory of elastic flow with isotropic hardening, the elastoplastic tangent is equal to where Ξ = C/ I + ∆λC : ∂ 2 σσ 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: 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.

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.

≈
At the same time, from the theory of elastic flow with isotropic hardening, the elastoplastic tangent is equal to = Ξ − Ξ: ⊗ : Ξ : Ξ: where Ξ = /( + Δλ : ∂ ). 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: where = Ξ: ⊗ : Ξ, and = : Ξ: (top index corresponds to each problem). Thus, by solving three problems with different boundary conditions (12), the soughtfor hardening parameter values can be obtained.

Research Object
A two-dimensional concrete structure with inclusions in the form of fibers is considered for different cases of fibers: 1. Evenly distributed and unidirectional; 2. Bidirectionally distributed; 3. Randomly distributed.
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.

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:

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: 1.
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 Crystals 2021, 11, 776 6 of 26 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. 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.
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  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.
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 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.
with full geometry (Figure 5a-c). This solution is considered as an exact solution. Th rors of the compared models are calculated referring to it. A coarse grid (Figure 5d quires a solution made using the effective elastic tensor. The models using the effe hardening factor and the concrete hardening parameter are compared. When using re mapping, the value of the plastic strain tensor for calculating the effective coefficie taken from the previous iteration. The computational mesh of the complete geometry for unidirectional, bidirecti and random fiber distribution and the coarse computational mesh, on which the solu is performed using the effective parameters, have 33,963 vertices, 78,361 cells, 44,814 tices, 100,063 cells, 47,822 vertices, 106,910 cells, and 2319 vertices, 4640 cells, respecti

Calculation of Effective Parameters
To calculate the effective parameters, for the unidirectional location of fibers, a cell is used. For the bidirectional location of fibers, we used a 2 × 2 cell. In the case w random distribution of fibers, the accuracy of the elasticity problem highly depend the representative volume size. It is thought that a 16-fiber case would be optimal. computational meshes for calculating the elastic parameters are shown in Figure 6a-f the computational mesh for calculating the hardening parameters is shown in Figure The computational mesh for calculating the elastic parameters for unidirectiona directional, and random distribution and the coarse computational mesh for calcula 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.

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.  After applying the algorithm, the effective parameters presented in Table 2 were tained. Table 2. Calculated effective elasticity and plasticity parameters for all cases: unidirectional, b rectional, and random distribution of steel fibers.     Table 1 for unidirectional fiber location case.  Table 1 for unidirectional fiber location case.     The distributions of the deformation error modulus |u err |, which can be expressed by fine mesh solution u f ine taken as exact and homogenized solution u hom as |u err | =  (a) (b) 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.
(a) (b) 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.  (a) (b) 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.
(a) (b) 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 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. 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.

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. 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.

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. (a) (b) 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.

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, 19, 20, 21, 22, respectively. 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.      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.

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. 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.

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 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.

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 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.

Figure 25.
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 random fiber distribution case.
The distributions of the deformation modulus error for test 3 for unidirectional, bidirectional, and random distribution of fibers are shown in Figures 26, 27, 28, 29, 30, respectively.     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.

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. 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.

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 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.

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, bidirectional, and random distribution of fibers are shown in Figures 34, 35, 36, 37, 38, respectively.      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 elasticplastic flow of the composite material. 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.

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.

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.