Evolution of the stiffness tetrad in fiber-reinforced materials under large plastic strain

The main aim of this work is to track the evolution of the stiffness tetrad during large plastic strain. Therefore, the framework of a general finite plasticity theory is developed. Some special cases are examined, and the case of a material plasticity theory is considered more closely. Its main feature is that the elasticity law changes during plastic deformations, for which we develop an approach. As sample materials, we use three types of fiber-reinforced composites. For numerical experiments and verification of the model’s predictions, finite element simulations of representative volume elements for uni-, bi- and tri-directional reinforced materials with periodic boundary conditions are used. From these, we extract the stiffness tetrads before and after large deformations of the material. We quantify the change of the stiffness tetrads due to the fiber reorientation. Finally, we propose an analytical evolution with three parameters that account reasonably well for the evolution of the stiffness tetrad.

In other cases, the assumption of a constant elasticity reference law is not applicable, for example, when the elastic directors deform with the material. The setting in which the directors of elasticity move with the material is referred to as material plasticity after [4,11] due to the fact that the elastic directors behave like material line elements. The most obvious example is fiber-reinforced materials. For the orientation distribution evolution of short fibers in viscous flow, theories have been developed, starting from the seminal work [15] to industrial application [21].
Mathematically, the difference to the usual large strain plasticity theory with a multiplicative decomposition of the deformation gradient is that the evolution of the stiffness tetrad and of the stress-free placement during plastic deformation is decoupled. To work out the details, we motivate and extend the usual equations by a second-order tensor which accounts for the evolution of the stiffness tetrad. We also compare this extension to the isomorphy concept after [2][3][4]10,23,28,29,31]. They have developed a finite plasticity theory which bases on the transformation of the elastic law during yielding. For many materials, it is known that the elastic behavior remains almost unchanged during large plastic deformations in the sense that in every elastic region after any plastic deformation the local elastic law is obtained from a push forward of the elastic reference law. Examples are crystal plasticity and isotropic von Mises plasticity. In both ways, we end up with the same equation for a general plasticity theory. We identify four special cases for this equation, one being the equation of our material plasticity theory.
Next, we like to compare to some other similar approaches from the literature. Boehler [5] proposed an invariant formulation of anisotropic constitutive equations. Within these equations, several scalar-valued functions appear which have to be identified in experimental investigations. In our work, we have to identify only three scalars. They can easily be determined using our RVE calculations which is rather simple and fast. Pastor et al. [24] use the formulation of Boehler [5] to define isotropic problems associated with originally orthotropic ones. Again, it is necessary to determine anisotropy coefficients and equations to define the association between the isotropic and the anisotropic problems. Menzel and Steinmann [22] started with the multiplicative elastoplasticity. Additionally, we also used the isomorphy concept as a second access to our equation. They introduce additional symmetric arguments into the scalar-valued isotropic tensor functions and develop a modular framework of the fundamental covariance relation. Again, a number of symmetric tensorial fields of second-order and scalar-valued internal variables are necessary which, in contrast to our approach, leads to much greater identification effort. Lu and Papadopoulos [19] also introduced a covariant formulation of anisotropic finite plasticity. Theoretical results are carried out, and structural tensors for characterizing the anisotropy are introduced. In contrast to our work, the evolution of the anisotropy during the deformation is not mentioned. In [16,17], the covariant formulation of Lu and Papadopoulos for describing the evolution of the anisotropy is mentioned. Tey uses a convected-type evolution equation for the structural tensors and an additional constitutive equation for the plastic spin. A close match-up with the experimental results of Kim and Yin [18] was reached. In [14], an elastoplastic constitutive model for orthotropic materials at large strain is presented. Within the theory, an evolution of the anisotropy axes due to rigid rotations is modeled. They also use the multiplicative decomposition of the deformation gradient and a quadratic yield function formulated in terms of the Mandel stress tensor.
In our work, we will develop the framework of a general finite plasticity. Some special cases are examined, and the case of a material plasticity theory is considered more closely. Its main feature is that the elasticity law changes during plastic deformations, for which we develop an approach.
To identify and compare the stiffness tetrads before and after large plastic deformations, a representative volume element (RVE) [13] with a fibrous microstructure is generated. Uni-, bi-and tri-directionally reinforced samples are considered using periodic boundary conditions. On the micro-scale, a standard isotropic elasticplastic material model without hardening is used.
After calculating the effective stiffnesses of the different material samples, we investigate their evolution during different deformations. It turns out that it is possible to denote the change of the stiffness with sufficient accuracy using one additional second-order tensor. After investigating the change in the stiffness tetrads, we finally propose an analytical evolution equation for this tensor.

Notation
A direct notation is preferred. Vectors are denoted as bold minuscules like v = v i e i . Summation over repeated indices is implied. e i are orthonormal base vectors for which e i · e j = δ i j holds, with the Kronecker symbol δ i j . The dyadic product and scalar contractions are denoted by with a dot being the usual scalar product between vectors. Second-order tensors are written as bold majuscules, like T = T i j e i ⊗ e j . The scalar contractions are carried out such that the positive definiteness is inherited to tensors, i.e., A : A = A i j A i j ≥ 0. When the product is clear from the context, as in a = Ab, the single scalar dot is occasionally omitted. In the case of multiple scalar contractions, as in T = K : E, we will always write the dots for better comprehensibility.
Normalized Voigt notation Fourth-order tensors play a central role in this work. They are denoted by blackboard bold letters, like the stiffness tetrad K. It has the subsymmetries and the major symmetry, which implies where the components are given w.r.t. the base tetrads e i ⊗ e j ⊗ e k ⊗ e l . A more convenient representation is obtained by introducing the symmetric base dyads can be arranged in a symmetric 6 × 6 matrix, for which the usual rules of matrix calculus hold since the basis is orthonormal. The normalization has been used implicitly firstly by [30], rediscovered by [20], and popularized by [32] in its non-normalized form, see also [6,9]. Other bases are not used in this work. Rayleigh product We further make use of the Rayleigh product, which is defined between a second-order tensor and a tensor of arbitrary order, It changes the basis while keeping the components, which can be used, e.g., to represent rotations. For secondorder tensors, it can be written as A * T = AT A T . It is associative in the first factor, i.e., ( AB) * K = A * (B * K), and linear in the second factor, i.e., A * ( Differences between tetrads Differences between stiffness tetrads are normalized w.r.t. K 1 , where K 0 is the initial stiffness and K 1 the stiffness after plastic deformation. Depending on the context, K 1 is the result of a numerical RVE experiment or a phenomenological model prediction, for which K 1 can be replaced by an expression with K 0 . To quantify a deviation, we take the Euclidean norm. The Euclidean norm of a fourth-order tensor is calculated by the square root of the fourfold contraction of the tensor with itself. For example, the normalized difference between a phenomenological model prediction and an RVE result is given by Another interesting result is the change of stiffness due to plastic deformations, i.e., the difference between initial and current stiffnesses.

Plasticity theory with the evolution of the stiffness tetrad
As explained in introduction, the aim of this work is to develop a phenomenological finite plasticity theory that describes the evolution of anisotropy. To explain this idea, a comparison with the known crystal plasticity theory is helpful. In the latter case, under plastic deformations, the crystal structure develops differently than the material (see Fig. 1a). In the theory of material plasticity, the opposite occurs. If there is a plastic deformation in the fiber and the matrix material, the fibers will deform together with the material and not differently. This means that the anisotropy induced by the fibers, in contrast to the crystal lattice, is firmly bonded with the material. The best example here is a fiber-reinforced material (see Fig. 1b).  For creating a phenomenological model for a material plasticity theory, we start with the kinematics. We define the reference placement, the stress-free placement and the current placement shown in Fig. 2. The multiplicative decomposition of the deformation gradient F (see, e.g., [1,7]) is Next, we choose a general elastic law in the stress-free placement with a material stress tensor mat T , the stiffness tetrad of the material as a tensor of fourth-order K and a strain measure. We use the Seth strain tensors with "m" being a real parameter and "e" marking the elastic part of a variable [26]. A general framework can simply be set up by postulating quadratic elastic strain energy in terms of the Seth strain tensors, We assumed small elastic deformations, such that it is sufficient to account only for the quadratic term in the strain energy. It would be possible to present our proposal with a general potential. A generalization could be obtained by a series expansion of the strain energy: One could then approach the change of all higher-order K tensors similar to our strategy for K (4th order) given below. It would be interesting to assess the quality of such an approach, but it this beyond the scope of our work. The elastic law is We write the elastic part C e as We now insert this equation into Eq. (8) and expand the I with F p to If we confine to the special case m = 2, it is possible to factor out F p , and thus, we can further simplify to The choice m=2 results in further simplifications. The stress tensor mat T becomes 2PK T SF in the stress-free placement. This also allows for a simple conversion to the Cauchy stresses T . The next step is to introduce our ansatz for capturing the evolution of the anisotropy axes during the plastic deformation. We know that the anisotropy axes transform by vectors. We will call the initial vector a 0 and the vector of the transformed axis a. Second-order tensors transform tangent vectors that represent line elements. We will name this tensor P K with Let us consider a transversal isotropic material. From representation theory, we can write the material stiffness tetrad K ti of this material as where "sym" is the symmetrization according to Eq. (1). Because we know this representation and there is only one vector a, we can transform the whole stiffness tetrad by transforming the vector a. In other cases, we either do not know the representation theorem or it changes during plastic deformation due to a change of the symmetry class (e.g., cubic symmetry). Because our theory has to be valid for all possible material symmetries, we have to make a simplification. In Eq. (16), we identify the summand with c 3 as the part capturing the main anisotropy features. The structure is similar to the Rayleigh product so we decide to directly apply the transformation P K to the whole stiffness tetrad K P K is defined in the stress-free placement and transforms the stiffness tetrad. It is not a change of placement. Furthermore, it is now possible to summarize F p and P K in the following equation and compare it with other theories. In the stress-free placement, the second Piola-Kirchhoff stresses are With the transformation between the stress-free placement and the reference placement using F p and its we can factor out F p and summarize with P K . We consider only isochoric plastic deformation and therefore J p = 1. Finally, we end up with the following equation for the second Piola-Kirchhoff stresses in the reference placement

Derivation of the theory with the isomorphy concept
We want to find another access to Eq. (20) and start with a physically linear elastic law according to St.Venant-Kirchhoff Denoting k 0 as a constant elastic reference law and k p as the current elastic law describing the material at any given time, the isomorphism P holds for the elastic law of an elasto-plastic material Substituting Eq. (21) into Eq. (22), one obtains the generalized form of a plasticity theory with isomorphic elastic ranges according to [4, p.285ff]. P is also called a plastic transformation and describes the change of the stress-free placement (C u0 ) and the stiffness tetrad (K 0 ) For our theory, the constant elastic reference law has to be replaced. A generalization of the isomorphy concept is to choose two different plastic transformations P K and P C , where P C transforms the stress-free configuration and P C P K transforms the elastic stiffness into the stress-free configuration: With the relation we obtain Eq. (20) for our general plasticity theory. In this equation, the isomorphy condition no longer holds because there are two different plastic transformations. With Orth + defining the special orthogonal group (rotational group with determinant +1) and P K ∈ Orth + , we are able to combine both transformations to one. This is always possible for crystals where P K ∈ Orth + always holds. We now want to examine this equation for possible special cases.

Special cases for the general equation
In this subsection, we use the special unimodular group "Unim + ," defining all tensors with determinant equal to +1. The special linear group "Inv + " defines all invertible tensors with positive determinant. By choosing P K = P C = I, it is clear that this will lead to an elasticity theory which can be used for elastic materials or the elastic range of a material (see Table 1, column 1) In the case that P K = I and P C ∈ Unim + , we will obtain the isomorphy concept which is, e.g., applicable as crystal plastictiy theory as described in [4] where only one plastic transformation exists (see Table 1, column 2) By choosing P K ∈ Inv + and P C = I, one would describe evolving material properties without a deformation of the material. This could be applicable to the description of some recrystallization processes where P K = Q (see Table 1, column 3) Q is an orthogonal tensor of second order and refers to a rotation. We refer to a recrystallization process on the micro-scale, where each material point has just one lattice orientation. This lattice orientation can change Evolving material properties Different evolution of material properties and stressfree placement Elastic materials Materials with crystal structure Fiber-reinforced materials by recrystallization, namely when grain boundaries move upon coarse grain annealing. The local orientation change can be expressed by a single orthogonal tensor. On the macro-scale, the orientation distribution evolves, which is, of course, not captured by a single second-order tensor. The last possible combination is to allow for a different evolution of P K and P C and thus a different evolution of material properties and the stress-free placement. In this case, one can describe, e.g., fiber-reinforced materials with the idea of a material plasticity theory explained in the following (see Table 1, column 4). The final equation for our theory is

Representative volume elements
In this section, we present the representative volume elements (RVEs) used in this article. We modeled unidirectional, bi-directional and tri-directional fiber-reinforced materials all having different material symmetries. These are purely academic examples with the aim of analyzing the development of the stiffness tetrads. Of course, real materials look different. To create the meshed RVE with the finite element software Abaqus, various individual steps are necessary. We used the scripting language Python which is able to access Abaqus functions and developed a script to generate the RVEs. All simulations are displacement controlled by prescribing the components of the effective displacement gradient H and using periodic boundary conditions as described in [12]. In all simulations, we used the 20 nodes quadratic elements C3D20. We also ensured that no large effective volume dilatations occur during the FE simulations by • prescribing only isochoric effective final deformations, • allowing for free lateral straining such that no volumetric straining is enforced, • prescribe a strain path without intermediate large volumetric strains that may occur when ramping up the effective displacement gradient linearly.
The RVE has two phases, namely the matrix and the fiber material. For both materials, we use an elasticplastic material model for large deformations using the isomorphy concept as described, e.g., in [4]. The isotropic elastic-plastic material has the material parameters Young's modulus E, the Poisson's ratio ν and the yield stress σ F . The three parameters have different values for the matrix and the fiber material. The model is implemented as UMAT subroutine into Abaqus. The elastic behavior is described by the St. Venant-Kirchhoff law. We have an isotropic yield condition following the J 2 theory by von Mises and a flow rule using the principle of maximum plastic dissipation. We do not incorporate compressible materials in our work. The main focus of our work lies in the anisotropy. Since volume changes are isotropic deformations, they do not contribute to the anisotropic changes. We also neglect viscosity, relaxation and the like. This is a strength of our method and a necessary simplification. If we mix too many effects, it will not be possible to identify individual contributions to the overall behavior. Here, we are concerned with the evolution of the elastic anisotropy due to a fiber reorientation due to plasticity. After having identified a reasonable approach and finding the parameters, it is possible to construct a more sophisticated model, for example, by adding a viscous contribution or hardening. But for identification purposes, we need to consider only the effect that we seek to model. The following sections show the different RVEs. The first RVE is supposed to have an initially effectively transversal isotropic material symmetry. Therefore, we choose an uni-directional reinforced material with circular fibers in a hexagonal arrangement in z direction. In Fig. 3, the unit cell is shown. The parts depicted in Fig. 3 are defined in dependency of the fiber volume fraction. Additionally, in Fig. 4 the meshed RVE is depicted in an undeformed and a deformed state.

Bi-directional reinforcement
The second RVE will have a tetragonal symmetry. For the sake of simplicity, the RVE for this material will have the shape of a cube. The cube is bi-directionally reinforced in the direction of x and y with fibers that have a square cross section. Figure 5 shows the unit cell. The generation of the meshed RVE is similar to the procedure of the uni-directionally reinforced RVE. In Fig. 6, the meshed RVE is depicted in an undeformed and deformed state.

Tri-directional reinforcement
The third RVE will be tri-directionally reinforced and exhibit initially a cubic symmetry. The RVE is nearly the same as for the bi-directional reinforcement, but now there are fibers three perpendicular directions. In Fig. 7, the unit cell is shown. In Fig. 8, the meshed RVE is depicted in an undeformed and a deformed situation.  In this section, we present a method to determine the effective stiffness of arbitrary inhomogeneous materials before and after large deformations. We will use the defined RVEs (see Sect. 3) to perform test deformations G E and use the resulting 2PK T stresses to calculate the stiffness tetrad K during the post processing. To track the evolution of the stiffness tetrads K, we have to determine the stiffness of the undeformed material (K 0 ) and the stiffness tetrad of a deformed material after a large plastic deformation (K 1 ). The effective stiffness tetrads will be calculated and compared in the reference placement, and therefore, we use the effective second Piola-Kirchhoff stresses 2PK T and the effective Green strain tensor G E. If the theory of material plasticity holds, there will be a second-order tensor P K which transforms the stiffness K 0 of the undeformed material to K 1 after a large plastic deformation.
For calculating the stiffness tetrads, we use the difference quotient out of six test calculations. Within the test calculations, six different small elastic test strains δ i with i=1,…,6 are applied to the presented materials (RVEs). With the definitions • H 0 : displacement gradient of the unloaded placement • H i : displacement gradients of the 6 different elastic deformations and the scheme in Fig. 9, we are able to calculate a stiffness tetrad by the equation The scheme is carried out at the (average) stress-free placement. To get to this placement, we unload componentwise and obtain the macroscopic P C . It turns out that due to the small elastic strains the difference between P C and the deformation gradient F −1 is very small. This simplifies greatly the empirical approach, since we can in essence avoid an independent evolution of P C but take P C = F −1 instead. The equation is implemented in a Mathematica script which comes with this article as supplementary material.

Test cases
We specified different material parameters for the matrix and the fiber material as shown in Table 2. These parameters were chosen as an academic example. Ten different test deformations were examined as depicted in Table 3. Whenever a strain component is not specified, it is implied to adjust freely such that the corresponding effective reaction stress is zero. In Eq. 31,  Table 3 Overview of the deviation ||K 1 − P C P K * K 0 || / ||K 1 || in % of the stiffness tetrads for the four suggested choices of P K ( P −1 C , I, P num 3 and P num 9 ) and the analytical P ana 3 for the three RVEs. The value of γ and is each 0.5. (Color figure online) we chose three elongation tests In this case, we always deform the material parallel or perpendicular to the fiber direction. With the free lateral straining, this corresponds to a uniaxial tensile state. Next, we choose six shear tests which lead to different They can lead to a transport of the fibers with and without an effect on the symmetry and to a change of the fiber direction. We will discuss this in the next subsection. The last test (33) combines one tensile and two shear tests and acts as the worst case

Results
We want to show some results for the uni-directional material. First, we perform a tensile test in fiber direction with 50% nominal strain and compare the values in the stiffness tetrads (given in GPa). We see in the case of uni-axial tension in Eq. (34) that there is no difference between K 0 and K 1 .
The second example is a shear test with H xy = 0.5 perpendicular to the fibers (Eq. (35)). We see that there is a difference between both stiffnesses, which implies that for certain deformations the stiffness tetrad evolves.
The reason is that we change the material symmetry during the deformation. The following three examples will show the deviation between the initial stiffness tetrad K 0 and the stiffness tetrad K 1 measured after a large deformation. We will investigate the three characteristic shear tests H xy , H xz and H zx . We see that the deviation of the stiffness tetrad K 1 from K 0 can be considerable. A shear parallel to the fiber normal plane with γ = 0.5 results in a deviation by approximately 35%. In the next section, we propose an analytical approach for predicting the evolution of K.
The first case is the shear test with H xy . Figure 10 shows that we have a transport of the fibers which leads to a change of the symmetry. The fiber direction remains the same. This type of deformation leads to a change of about 8% between the stiffnesses.
In Fig. 11, we see the shear test H xz = 0.5 and no transport of the fibers. Instead, the fiber direction changes which leads to a change of the symmetry. We find the largest change in the stiffness tetrad as shown in Eq. (37).
The shear test H zx = 0.5 is depicted in Fig. 12 and shows that we have a transport of fibers in the fiber direction. The fiber direction remains the same so this motion has nearly no effect on the symmetry

Analytical evolution equation
In Sect. 4, we presented a method to determine the stiffness tetrads before and after large plastic strains. The results show that the stiffness tetrad changes significantly, and that this change cannot be predicted by the plastic transformation P C of the isomorphy concept. Therefore, we decide to use the theory of material plasticity (see Sect. 2) using two different plastic transformations P C and P K . When we think of the form of the evolution laws for P K and P C , we need to distinguish the scales. On the micro-scale in the RVE, we locally use the von Mises flow rule, which is thermodynamically consistent. On the macro-scale, the empirical evolution of P K alone that we propose is not enough to assess the thermodynamic consistency. Due to P C = F −1 , i.e. F e = I, we have a rigid plastic process in which, somewhat paradoxically, we track the evolution of the stiffness indirectly by tracking the evolution of P K . In the present form, the approach is only applicable when the strain process is given. This may be enough for estimating the stiffness evolution in a forming process. To use the approach in a boundary value problem, we have to drop the simplification P C = F −1 to be able to calculate stresses and search the equilibrium. In that case, we need to complete the set of equations by adding a flow rule for the evolution of P C , for example by flow in direction of maximum plastic dissipation. Then thermodynamic consistency may be checked. On a side note, we believe that the approach is thermodynamically safe for reasonable flow rules. In any case, errors due to the potential violation of thermodynamic consistency are probably much smaller than errors due to the overall simplicity of the approach.

Evolution of the transformation P C
Following our assumptions, we have only small elastic strains. This directly leads us to the major simplification that This simplification has already been underpinned by the results in Sect. 4. The difference between P C and F −1 is very small and a result of the difference of order of magnitude between the yield stress σ F and the Young's modulus E in the two phases.

Finding a suitable transformation
The aim is to finally develop an evolution equation for the transformation P K of the stiffness tetrad K. To find this evolution equation, we want to study different possibilities. The first and one of the easiest approaches will be to choose P K = P −1 C . This would mean that the stiffness tetrad does not transform during a plastic transformation in contrast to the unloaded placement which will transform accordingly to P C . The second choice will be P K = I. This is the known isomorphy case and we already know that this approach will fail. For small plastic transformations, we can expect only a slight deviation while for large plastic deformations it will increase. The third possibility, we want to study, is to numerically fit a second-order tensor P modifying the known transformation P C . We will use two different sets of nonzero components. The numerically best fit is calculated with a minimization procedure in Mathematica. In the first case, we choose to take all possible nine components of a second-order tensor and call it P K = P num 9 with P num 9 = ⎡ ⎣ P x x P xy P xz P yx P yy P yz P zx P zy P zz If the scalar contraction of P C with a second-order tensor is able to track the evolution of the stiffness tetrad, we will get the exact results for the numerical fit of this choice. The second choice is to consider only the main tension Hxx=0.5, shear Hxy=0.5, Hxz=0.5 shear Hxy = 0.5 tension Hxx = 0.5 Fig. 13 Deviation of the measured stiffness tetrad from the respective modeling approach for the tri-directional reinforcement. The choice P K = P −1 C is clearly the worst, followed by elastic isomorphy P K = I. With P num 3 having three degrees of freedom, a significant improvement is reached. The improvement using P num 9 with 9 free values is not much better. The analytical approach P ana 3 is almost as good as the corresponding numerical fit P num 3 with 3 parameters. (Color figure online) diagonal components and call this transformation P K = P num This means we limit the field of application to materials having their anisotropic axes perpendicular to each other, where the directions of anisotropy are along e i . The main point of this idea is that P K transforms the original stiffness tetrad K 0 in the reference placement. Only afterwards, the transformation P C is applied.
To compare the different choices, we choose the tri-directionally reinforced material as example material. The material properties are the same as in the foregoing sections. We perform a shear test, a tensile test and the combined test of two shears and one tension and compare the deviation of the initial stiffness K 0 and the stiffness K 1 after the deformation (see Sect. 4). The deviation of the measured stiffness tetrad from the respective modeling approach for the tri-directional reinforcement is depicted in Fig. 13. The choice P K = P −1 C is clearly the worst. There is no agreement if we do not transform the initial stiffness tetrad after such large deformations. With the approach of elastic isomorphy P K = I, we find still a large deviation for the shear test and the combined test. In this case, the anisotropy of the stiffness tetrad is changed during the deformation. During the tensile test, the choice of P K = I works because during this deformation, the symmetry of the material is not changed. Using the transformation P K = P K3 with the three degrees of freedom, a significant improvement is reached. The improvement using P num 9 is not much better. The main results of this section are: Firstly, it turns out that it is possible to approximate the evolution of the stiffness tetrad with sufficient accuracy using one second-order tensor. Secondly, we find the best compromise between accuracy and effort with P K = P num 3 .

Evolution of the transformation P K
To derive an analytical evolution equation forṖ There is no direct physical interpretation for the pull-back of the velocity gradient, just like there is no direct physical interpretation for the reference placement or any other quantity that is pulled back. Nevertheless, they are mathematical auxiliary quantities that may be used. The main purpose of such pulled back material quantities is to formulate constitutive laws in accordance with the principles of material modeling, specifically the principle of invariance under superimposed rigid body motions. For this reason, the St. Venant-Kirchhoff law is used, with K 0 defined in the reference placement. Consequently, when modifying the elastic law, we need to do it in the reference placement, but of course, the actual deformation determines its evolution. Hence, we need to realize some kind of pull-back of a kinematic variable to model the evolution of K 0 . We could also have used a push-forward of the elasticity law where then the evolution of F * K 0 is tracked, but then we would mix the change of K due to rigid rotations with the change due to the deformation. It is precisely for this separation of effects that we describe the evolution of K 0 in the reference placement.
Since we want to make the approach as simple as possible, it should be linear. In Sect. 4.3, we have seen that shear has the greatest influence on the anisotropy change. Looking at a single fiber, there is always transverse isotropy. Hence, the shear directions in the cross-sectional plane are equal; we determine the shear rate by the theorem of Pythagoras to get the amount of the shear rate. To derive the evolution of the three main diagonal entries ofṖ we take the Euclidean norm of the corresponding components of L ref following Eq. (45) We have to determine the coefficients once for the chosen material. Therefore, we use RVE calculations for the three characteristic shear tests H yx (k 1 ), H xy (k 2 ) and H xz (k 3 ) and take the i'th entry in the numerically best P num 3 for k i . We verified the analytical evolution equation of P K by using the new strain path of the mixed test H xy , H xz and H x x . Afterwards, we used these constants to calculate all other tests. In Fig. 13, we finally show the results of the analytically calculated P K = P ana 3 . It turns out that the analytical solution is very close to the numerical best P num 3 . This is a highly satisfying result because the assumptions and simplifications we made were extensive. First, we assumed that a second-order tensor connected with the Rayleigh product to the fourth-order stiffness tetrad is enough to capture the evolution of the stiffness tetrad during a symmetry changing deformation. Secondly, we reduced the number of free components within the second-order tensor from nine to three in order to be able to find an analytical evolution equation which is easy to access. Our approach clearly gives better results than the known isomorphy concept and the multiplicative decomposition where the material stiffness transforms accordingly to the stress-free placement. In Table 3, we provide the results of all three different materials and ten different tests.
Finally, we like to represent and discuss the evolution of the stiffness tetrad during the transformation with the chosen second-order tensor P ana 3 . For this purpose, we use the shear test H xz on the uni-directional reinforced material as depicted in Fig. 14.
We will analyze the 21 different components of the stiffness tetrad itself and additionally the six Eigenvalues of the stiffness tetrad. In Fig. 15, we see the change of the components during the shear test. Additionally, we show the index arrangement of the Voigt notation for a stiffness tetrad K, the absolute values of K 0 and K 1 in Fig. 16.
Because the shear test is in the 1-3 planes, we see the largest change in the components "xxxx" and "zzzz" marked in dark red, in the components containing two times "x" and "z," namely "xxzz" and "xzxz" marked in red and in the two components "xxxz" and "zzxz" marked in light red. All other components containing at least one time the second index are marked in blue. In dark blue, we marked the components with a large change of nearly 100%. These components were zero in the initial stiffness K 0 and changed to some big value in K 1 . The components are "yyxz" and "xyyz.' The components with a medium change are "xxyy." "yyzz," "xyxy" and "yzyz" are marked in blue. The component "yyyy" has an initial value but zero change and is marked in light blue. All values with a zero initial value and zero change are marked in black and are "xxxy," "yyxy," "zzxy," "xyxz," "xxyz," "yyyz," "zzyz" and "xzyz." In Fig. 17, we additionally investigated the evolution of the Eigenvalues during the deformation. We sorted them in descending order from the largest one (dark red) to the smallest one (dark blue). Because the Eigenvalues are free of any rotation, we clearly see that the stiffness tetrad evolves and the symmetry changes.

Summary
Finally, we like to summarize the article. In Sect. 1, we give an overview of similar works to help you better understand our article. We introduce the topic of our work and explain the structure of the following article.
In Sect. 2, we postulated a general equation for different plasticity theories. We explained two different ways to access this general equation. First, we start with the multiplicative decomposition of the deformation gradient. The second access to our theory is given by the isomorphy concept. We formulated four special cases for the obtained general equation and summarized them. One of the special cases is the so-called material plasticity theory. Here, we allowed a different evolution of the material symmetry within the stiffness tetrad and of the stress-free placement. We were also able to account for an evolution of the stiffness tetrad. The special case of fibers moving as material line elements is examined more closely.
We set up a numerical laboratory to investigate the material plasticity theory in Sect.
3. An elastic-plastic material model for large deformations using the isomorphy concept on the micro-scale was applied for the RVE calculations. We defined different material parameters within this model for the matrix and the fiber phase. For these calculations, we used the finite element program Abaqus with its scripting language Python. The computer algebra system Mathematica was used for the post-processing. Using three different RVEs for modeling an uni-directional, a bi-directional and a tri-directional reinforced material with different material symmetries, we were able to determine the effective stiffness of these materials in Sect. 4. Therefore, we implemented and validated a method to measure material stiffness tetrads using the difference quotient. The results lead us to our material plasticity theory as discussed in Sect. 2. We need a second-order tensor P K which transforms between K 1 and K 0 . It will predict the evolution of the stiffness tetrad during a plastic deformation of an anisotropic material which changes the material symmetry.
In Sect. 5, we found that it is possible to predict the evolution of the stiffness tetrad during large deformations using one second-order tensor P K within the equation of our material plasticity theory: We ascertained an analytical evolution equation to predict this tensor P K with only three free variables. This is valid for a class of fiber-reinforced materials having fiber angles of 90°with a sufficient accuracy. Finally, we like to give an outlook. With our theory, we are able to predict the evolution of the stiffness tetrad during large plastic deformations. For this purpose, we finally used a second order tensor with three free variables which is, for us, the best compromise between accuracy and applicability. For more accurate results or materials with a micro-structure different to our 90°scheme, one could use more free variables. One can also investigate the third special case within our general equation for plasticity theories. Here, it is possible to describe evolving material properties.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.