Linearization of Composite Material Damage Model Results and Its Impact on the Subsequent Stress–Strain Analysis

To solve problems in the field of mechanical engineering efficiently, individual numerical procedures must be developed, and solvers must be adapted. This study applies the results of a carbon-fibre reinforced polymer (CFRP) analysis along with the nonlinear finite element damage (FE) method to the translation of a linear solver. The analyzed tensile test sample is modelled using the ply-by-ply method. To describe the nonlinear post-damage behavior of the material, the Hashin model is used. To validate the transformation, an analysis and comparison of the damage results of the linearized and nonlinear model is carried out. Job linearization was performed by collecting elements into groups based on their level of damage and pairing them with unique material cards. Potentially suitable mathematical functions are tested for the grouping and consolidation of the elements. The results show that the agreement of some presented methods depends on the damage level. The influence of the selected statistical functions on the result is shown here. The optimal solution is demonstrated, and the most efficient method of linearization is presented. The main motivation behind this work is that the problem has not been discussed in the literature and that there is currently no commercial software translator that provides the transference of models between solvers.


Introduction
In the field of mechanical engineering, there are several problems that need to be solved by numerical methods. To increase their efficiency, individual procedures have been developed for individual tasks and solvers have been adapted. The most common variants of the solution include a linear or nonlinear variant of the calculation following an explicit or implicit scheme [1,2]. Individual calculation schemes differ in mathematical formulation, stability and speed of computation. Linear solvers are based on the assumption that the force in the model increases with increasing deformation during the calculation [1,2]. If this condition is violated, the task becomes poorly conditioned and the calculation cannot be successfully performed [3]. However, the phenomenon of force drop in the model occurs in all tasks that are focused on the study of the state of structures in a broken state [4][5][6][7] beyond the strength reached before total structural collapse [8][9][10][11][12]. This can occur, for example, due to a collision between fast-moving objects [13]. For these reasons, violation tasks are calculated with the help of solvers following the so-called explicit scheme. They do not require the condition of increasing forces in the model; they are stable but only suitable for simulations of fast or quasi-static processes. Because the tasks for individual solvers are formulated differently, it is necessary to ensure appropriate continuity between them. Thus, the order of calculations, which are based on different calculation schemes, cannot be created completely arbitrarily. A good example is the sequence of analyses to determine the stiffness of a composite part before and after failure [14], which is the focus of this article, because changes in stiffness can be based on a structural condition monitoring (SHM) system that can determine the damage location of a local stiffness change. First, a linear approach following the implicit scheme is used, then a nonlinear analysis with violation is performed using the explicit solver, and then again, a linear analysis with the implicit scheme follows. As an example, in this work, a method of transferring information from the results of the explicit solver to the model of the implicit solver Abaqus 6.14 is presented. The linearization method is used to transfer information between individual models. The main motivation of this study is that there is no communication between solvers and there are not enough resources to address a similar issue in such a direct way, as addressed in this work. Furthermore, a commercial software translator providing the translation of models does not currently exist. This paper focuses mainly on composite material failure, such as has been referred to in previous studies [15][16][17][18][19][20]. In general, the linearization (model conversion) procedure can be applied to all materials and implicit solvers. The aim of the article is to find a unique and highly efficient linearization method, which is theoretically verified.

Workflow Overview
The workflow diagram is shown in Figure 1. The experimental simulation following the standardized ASTM D3039 test was chosen for demonstration [20]. The problem contains a material description at the level of Hooke's law [1,3] and follows Hashin's progressive damage model [1,17,19]. The results are subject to translation, which mainly concerns the material description at the level of individual elements of the model. According to the degree of damage, the damaged elements are consolidated into groups with help of mathematical functions and paired with unique material cards, which are calculated for each group. Furthermore, an analysis and comparison of linear and nonlinear models is performed. The success of consolidation functions is compared in a demonstration task.
Polymers 2022, 14, x FOR PEER REVIEW 2 of 1 condition monitoring (SHM) system that can determine the damage location of a loca stiffness change. First, a linear approach following the implicit scheme is used, then a non linear analysis with violation is performed using the explicit solver, and then again, a lin ear analysis with the implicit scheme follows. As an example, in this work, a method o transferring information from the results of the explicit solver to the model of the implici solver Abaqus 6.14 is presented. The linearization method is used to transfer information between individual models. The main motivation of this study is that there is no commu nication between solvers and there are not enough resources to address a similar issue in such a direct way, as addressed in this work. Furthermore, a commercial software trans lator providing the translation of models does not currently exist. This paper focuse mainly on composite material failure, such as has been referred to in previous studies [15 20]. In general, the linearization (model conversion) procedure can be applied to all mate rials and implicit solvers. The aim of the article is to find a unique and highly efficien linearization method, which is theoretically verified.

Workflow Overview
The workflow diagram is shown in Figure 1. The experimental simulation followin the standardized ASTM D3039 test was chosen for demonstration [20]. The problem con tains a material description at the level of Hooke's law [1,3] and follows Hashin's progres sive damage model [1,17,19]. The results are subject to translation, which mainly concern the material description at the level of individual elements of the model. According to th degree of damage, the damaged elements are consolidated into groups with help of mathe matical functions and paired with unique material cards, which are calculated for each group. Furthermore, an analysis and comparison of linear and nonlinear models is per formed. The success of consolidation functions is compared in a demonstration task. In principle, the individual formulations differ in one specific way in terms of mate rial description. While the nonlinear violation model contains mathematical function with violation parameters, in the linear model, nonlinearities are replaced by a large In principle, the individual formulations differ in one specific way in terms of material description. While the nonlinear violation model contains mathematical functions with violation parameters, in the linear model, nonlinearities are replaced by a larger number of material cards with a linear response. The consolidation functions listed below provide the matching key between the material cards and the violation parameters.

Basic Calculus
In this paper, the components of the stress and strain tensor in the form shown by Equation (1) are considered for each finite element [1]. This corresponds to the formulation for tasks in the field of composite materials, solved with the help of the Abaqus system and with the use of shell-type elements. Components of tensors that are perpendicular to the plane of the element are not considered, so each contains only five components instead of the usual six.
The linear formulation of the stress and strain analysis problem is shown by Equation (2) [1,18,19], where [C 0 ] in Equation (3) represents the initial stiffness matrix of the composite material without including the possible failure effect.
The stress calculation of the damaged composite layer in simplified form is also performed using the Equations (6) and (7) [1], where C d represents a stiffness matrix with damage rates d 1 , d 2 , d 12 , d 13 and d 23 , similar to Equation (4); however, this formula does not include the influence of Poisson's ratio decrement, such as in (7) [9].
The damage rates in the main directions of anisotropy d 1 , d 2 and their relation to the effective stress show the Equations (9) and (10) [1,8].
Equation (11) shows that the stresses calculated using the linear and nonlinear formulation are equal up to the ultimate strength of the material.

Graphic Interpretation of Hashin's Violation
A graphical interpretation of the stress response of the composite material to deformation in the main direction of anisotropy is shown in Figure 2 [1,6]. The diagram does not express all the key points of Hashin's theory, instead showing only those that are essential in this paper and for the solver used [1].
The degree of shear damage can be expressed by the Equation (8) [1], where d , d , d and d represent the rates of damage in tension and pressure.
The damage rates in the main directions of anisotropy d , d and their relation to the effective stress show the Equations (9) and (10) [1,8].
Equation (11) shows that the stresses calculated using the linear and nonlinear formulation are equal up to the ultimate strength of the material.

Graphic Interpretation of Hashin's Violation
A graphical interpretation of the stress response of the composite material to deformation in the main direction of anisotropy is shown in Figure 2 [1,6]. The diagram does not express all the key points of Hashin's theory, instead showing only those that are essential in this paper and for the solver used [1]. Before failure, the linear solver follows the blue line without stiffness changes with material stiffness E , which is the same in tension and pressure. The nonlinear solver proceeds when breaking a composite material according to Hashin's theory, which is based Before failure, the linear solver follows the blue line without stiffness changes with material stiffness E 0 i , which is the same in tension and pressure. The nonlinear solver proceeds when breaking a composite material according to Hashin's theory, which is based on the assumption that the breaking can take place independently of the meaning, direction and plane of anisotropy. Once the material strength limit is reached, the material stiffness is gradually lost as a result of the failure. The accumulation of a larger number of failures then leads to a depletion of "damage energy" and a complete loss of rigidity. Thus, reformulating the problem for the linear solver has obvious pitfalls. There may be a situation where a total failure occurs, or the material is damaged in tension (d + i ) and pressure (d − i ) in different ways. This is reflected in all elastic constants. These can then be unique for each element or their group, depending on the degree of damage. The linear solver in post-damage analyses uses only one stiffness value and proceeds according to the curve E d i . In the presented method, we assume that we can categorize groups of elements into independent groups according to the degree of damage. Unlike the nonlinearity of elastic properties, the uniqueness of material properties is not an obstacle for a linear solver. When transforming a problem with a violation, it is therefore necessary to reformulate the results, especially in terms of material properties.

Consolidation Function
The damage of each element is characterized by four state variables: d + 1 , d − 1 , d + 2 and d − 2 , which represent its size, direction and meaning. The consolidation function has two roles. First, it must facilitate the association of state variables and the categorization of elements according to the degree of damage. The main goal of this work was to find and verify suitable consolidation functions F1, F2 and F12, with which it is possible to effectively perform parameter reduction. Functions were chosen and sorted with respect to the complexity of their formulation and the truth of the response. The state variables and their consolidation scheme are shown in Table 1.
A total of five functions were tested, which are listed in Table 2. These are basic statistical functions or functions that can be easily derived from them. Table 2. List of damage consolidation functions.

Method Consolidate Condition
Arithmetic mean The arithmetic mean was calculated using the Equation (12).
The weighted average was calculated using Equation (13), where d i represents the standard calculated according to (14). This is used in the calculation of the weights w + i a w − i in Equation (15).
The consolidation of shear damage is found using the function F12, which takes the form of form (8) for all methods.

Material Card Generation
This research shows that perhaps the most accurate way to translate the results of a nonlinear problem into a linear problem is to assign a unique material card to each element. The number of elements in the model then corresponds to the number of material cards. However, this approach places extremely high demands on the preprocessor, making the task confusing.
Therefore, the material cards were generated according to the diagram shown in Figure 3, which represents different variations of damage with a 10% step for each meaning and direction. The variation in damage in one direction can then be represented by an array of 11 rows and 11 columns. However, given the stability conditions (Table 3), a real material card cannot exist for zero 100% damage. At the extreme intervals, in the context of stability conditions [1,18], we must consider a maximum damage value of 99.99%. As the results show, this consideration is sufficient for demonstration purposes. The field then contains exactly 11 × 11 items, which represents a total of 121 material cards.
for i= 1, 2 The consolidation of shear damage is found using the function F12, which takes the form of form (8) for all methods.

Material Card Generation
This research shows that perhaps the most accurate way to translate the results of a nonlinear problem into a linear problem is to assign a unique material card to each element. The number of elements in the model then corresponds to the number of material cards. However, this approach places extremely high demands on the preprocessor, making the task confusing.
Therefore, the material cards were generated according to the diagram shown in Figure 3, which represents different variations of damage with a 10% step for each meaning and direction. The variation in damage in one direction can then be represented by an array of 11 rows and 11 columns. However, given the stability conditions (Table 3), a real material card cannot exist for zero 100% damage. At the extreme intervals, in the context of stability conditions [1,18], we must consider a maximum damage value of 99.99%. As the results show, this consideration is sufficient for demonstration purposes. The field then contains exactly 11 × 11 items, which represents a total of 121 material cards.  Table 4. Other shear modules can be calculated in an analogous way. The derivation of Equation (16) is given in the literature [18,19], where its validity is considered only within the material without damage. In this contribution, the validity of the relationship is extended to all material cards, including material cards with violations.   Assuming the validity of the conditions given for Equations (1)-(11), the elastic properties of the individual material cards can be calculated using (16)-(22), for which the properties are listed in Table 4. Other shear modules can be calculated in an analogous way. The derivation of Equation (16) is given in the literature [18,19], where its validity is considered only within the material without damage. In this contribution, the validity of the relationship is extended to all material cards, including material cards with violations. , if σ 11 = 0, σ 22 = 0, σ 12 = 0 (20)

Material Properties
For test tasks, a material with elastic parameters was used, which are listed in Table 5. Strength parameters respecting the rules of Hashin's violation in the Abaqus system are shown in Table 6. Parameters describing the progressive development of the violation are given in Table 7 [14].   Figure 3.

Damage Variation Matrix and Response of Consolidation Functions
The graphical display illustrates the consequence of the relationship of the selection group of elements at the level of one direction of damage. While Figure 4 presents the consolidation of the damage parameter in the direction of the main diagonal, Figure 5 shows the direction of the side diagonal. Figure 6 describes the variation of damage in the top row direction, and Figure 7 displays the direction of the left column of the variation scheme. It is obvious that the damage values can take values that are only in the interval from 0 to 1.  The graphical display illustrates the consequence of the relationship of the selection group of elements at the level of one direction of damage. While Figure 4 presents the consolidation of the damage parameter in the direction of the main diagonal, Figure 5 shows the direction of the side diagonal. Figure 6 describes the variation of damage in the top row direction, and Figure 7 displays the direction of the left column of the variation scheme. It is obvious that the damage values can take values that are only in the interval from 0 to 1.   The graphical display illustrates the consequence of the relationship of the selection group of elements at the level of one direction of damage. While Figure 4 presents the consolidation of the damage parameter in the direction of the main diagonal, Figure 5 shows the direction of the side diagonal. Figure 6 describes the variation of damage in the top row direction, and Figure 7 displays the direction of the left column of the variation scheme. It is obvious that the damage values can take values that are only in the interval from 0 to 1.

Testing Examples
The simulation of the laminate sample tensile test was carried out according to ASTM D3039. The sample dimensions of the experiment diagram are shown in Figure 8. At point A, the sample was completely fixed; at point B, a kinematic excitation was introduced, which follows the course shown in Figure 9.

Testing Examples
The simulation of the laminate sample tensile test was carried out according to ASTM D3039. The sample dimensions of the experiment diagram are shown in Figure 8. At point A, the sample was completely fixed; at point B, a kinematic excitation was introduced, which follows the course shown in Figure 9. The demonstration task was modelled in all cases with the help of S4R shell elements with a Lamina material formulation and a mesh size of 1 mm. The sample was modelled using the ply-by-ply method. In the case of the explicit solver, the layers were connected by a cohesive contact without failure; in the case of a linear solution, it was a linear contact. The demonstration task was modelled in all cases with the help of S4R shell elements with a Lamina material formulation and a mesh size of 1 mm. The sample was modelled using the ply-by-ply method. In the case of the explicit solver, the layers were connected by a cohesive contact without failure; in the case of a linear solution, it was a linear contact. In the first cycle, the sample was loaded in the linear region and the initial stiffness was calculated. The second cycle was designed to break the sample. This part of the task was solved only using Abaqus/Explicit. The third cycle was again simulated using both Load amplitude (mm) Step time (s) Test load cycle 1st cycle 2nd cycle 3rd cycle Figure 9. Point B motion history.
The demonstration task was modelled in all cases with the help of S4R shell elements with a Lamina material formulation and a mesh size of 1 mm. The sample was modelled using the ply-by-ply method. In the case of the explicit solver, the layers were connected by a cohesive contact without failure; in the case of a linear solution, it was a linear contact.
In the first cycle, the sample was loaded in the linear region and the initial stiffness was calculated. The second cycle was designed to break the sample. This part of the task was solved only using Abaqus/Explicit. The third cycle was again simulated using both solvers. Table 8 shows the laminate compositions on which the process was verified.
Although it is only a simulation of the tensile test, due to the contraction of differently oriented undamaged layers of the laminate, local significant pressure damage also occurs. The occurrence of this well-known accompanying phenomenon is checked, because otherwise the validation would only be partially valid. This is also the main reason as to why the method is tested on different layouts.

Results and Discussion
In Figure 10, the initial stiffness of samples determined by calculation is shown. The result is the same stiffness for calculations obtained from linear and explicit solvers. As shown in Figure 11, the results of the percentage decrease in the stiffness of the tested samples with respect to the reference stiffness.
The agreement of the residual stiffness results, compared with a reference value, is shown in Figure 12. The deviance from the reference solution is referenced in Table 9. It represents a comparison between the initial stiffness calculated with a linear solver and residual stiffness calculated by an explicit solver. It is visible that the decrease in stiffness depends on the composite layup, while the applied load cycle is the same. This was the aim, because the influence of fibre orientation has been considered. In this case, the result of the nonlinear simulation of cycle 1 was considered to be the reference value.

Results and Discussion
In Figure 10, the initial stiffness of samples determined by calculation is shown. The result is the same stiffness for calculations obtained from linear and explicit solvers. As shown in Figure 11, the results of the percentage decrease in the stiffness of the tested samples with respect to the reference stiffness.
The agreement of the residual stiffness results, compared with a reference value, is shown in Figure 12. The deviance from the reference solution is referenced in Table 9. It represents a comparison between the initial stiffness calculated with a linear solver and residual stiffness calculated by an explicit solver. It is visible that the decrease in stiffness depends on the composite layup, while the applied load cycle is the same. This was the aim, because the influence of fibre orientation has been considered. In this case, the result of the nonlinear simulation of cycle 1 was considered to be the reference value.
In Figures 13-16, the groups of undamaged laminate elements (Layup 4), created based on the use of the presented consolidation functions, are shown.
It is important to note that where there was no ( Figure 13) or complete damage, all consolidating functions except the arithmetic mean results were the same. This is consistent with Figures 4-7. Thus, it means that Layup 1 and Layup 6 and 7 should theoretically have the same results; however, this did not occur. This was caused by the results of the linearization of several elements which are not completely damaged (though may be partially damaged). This error cumulates and because the stiffness depends on the spatial distribution of the element, the deviation from the true value is very high. This small difference is shown in Figures 14c,d and 15c,d. It is visible that the number of elements is nearly the same, but only the cracks in Figures 14d-16d are correct and acceptable.    In Figures 13-16, the groups of undamaged laminate elements (Layup 4), created based on the use of the presented consolidation functions, are shown.     It is important to note that where there was no ( Figure 13) or complete damage, all consolidating functions except the arithmetic mean results were the same. This is consistent with Figures 4-7. Thus, it means that Layup 1 and Layup 6 and 7 should theoretically have the same results; however, this did not occur. This was caused by the results of the linearization of several elements which are not completely damaged (though may be partially damaged). This error cumulates and because the stiffness depends on the spatial distribution of the element, the deviation from the true value is very high. This small difference is shown in Figure 14c,d and Figure 15c,d. It is visible that the number of elements is nearly the same, but only the cracks in Figures 14d, 15d and 16d are correct and acceptable.

Conclusions
It has been shown that in the case of small-scale damage, the majority of the residual stiffness is formed by the undamaged layer; therefore, the agreement of the results is still very good. However, in cases of significant damage, the tested statistical functions fail. In all cases, only the product consolidation function can work properly. It was also found that the reduction in the amount of information caused by the idealization of material cards does not have a significant effect on the result, unlike the choice of the selection function. The pairing of material cards with selected groups of elements, which was carried out in the way described, is unique, highly efficient and accurate. The calculation was stable in all cases.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interests.

(MPa)
Compressive elastic of damaged material in the direction 2 d 1 (1) Consolidated damage factor in the direction 1 d 2 (1) Consolidated damage factor in the direction 2 d 12 (1) Consolidated damage factor in plane 12 d s (1) Consolidated shear damage factor d + 1 (1) Tensile damage factor in the direction 1 d − 1 (1) Compressive damage factor in the direction 1 d + 2 (1) Tensile damage factor in the direction 2 d − 2 (1) Compressive damage factor in the direction 2 Ultimate compressive strain in direction 2 U + 1 (mJ/mm 2 ) Tensile fracture energy in direction 1 U − 1 (mJ/mm 2 ) Compressive fracture energy in direction 1 U + 2 (mJ/mm 2 ) Tensile fracture energy in direction 2 U − 2 (mJ/mm 2 ) Compressive fracture energy in direction 2 F1 Consolidation damage function in direction 1 F2 Consolidation damage function in direction 2 F12 Consolidation damage function in plane 12