Modeling of Compression Test of Natural Fiber Composite Sections

In recent years, there has been an increasing interest in the composite materials reinforced with natural fibers. Due to the easy and cheap methods of obtaining raw materials, the possibility of recycling, biodegradability, production and processing safe for health, such materials can be a good alternative to the composite materials reinforced with glass or carbon fibers. However, due to the lower mechanical properties of natural composites, their use as construction materials is still limited. Nevertheless, natural fiber composites have the characteristics that can be used in structural applications as long as the mechanical behavior is well understood, reliable and predictable. The paper presents the results of numerical calculations of the compression process of a composite reinforced with a fabric made of flax and jute fibers (trade names: Biotex Flax 400g/cm2 and Biotex Jute 400 g/cm2) on a basis of Kinetix R240 epoxy resin. The data necessary for the numerical analysis were calculated in the Digimat software using the Double Inclusion micromechanical model, while the simulations of compression of the details were carried out in the Ansys software. Sections with different number of layers were tested. The results were compared with the experiment. The buckling forces obtained in the numerical analysis are comparable to the experimental results. Two types of C – section buckling modes were obtained and they consisted of two or three half-waves of buckling.


INTRODUCTION
Nowadays, there is a noticeable increase in the use of plant-origin fibers as fillers in the polymer matrix. The main fillers of plant origin used in the form of short and continuous fiber are: hemp, jute, flax, sisal, pineapple, etc. These materials are characterized by a relatively low density while maintaining high strength and stiffness. It should be noted that these types of fillers are fully renewable and their production requires little energy expenditure. On the other hand, their disposal does not generate toxic fumes/waste. These fibers are also characterized by lower durability than in the case of synthetic composites, higher water absorption, lower degradation temperature and high variability of properties, which directly affects the high variability of the properties of composites filled with this type of fibers. The positive aspects of the above-mentioned fibers significantly outweigh the above-mentioned disadvantages; hence, in the current world trend to reduce the storage of synthetic waste, these materials will be still improved and used to an increasing extent [1÷3]. In most cases, the mechanical properties of composites are improved by adding fibers to the polymer matrix, because the fibers have much higher strength and stiffness than the polymer matrix [4,5]. On the other hand, a significant problem in the use of natural fibers is the fact that the water absorption increases along with the proportion of fibers in the polymer matrix, which may consequently lead to the destruction of the composite [6]. However, it should be noted that the properties of fibers of plant origin are characterized by a large discrepancy [7,8]. For example, in the literature it was mentioned that the Young's modulus of cotton fibers ranges from about 5.5 to about 12.6 MPa, while the tensile strength of ramie fibers can range from about 400 MPa to about 938 MPa [9]. The variability of the properties of a particular type of fiber may result from variable geometry of the fibers, structural and material properties. Interestingly, the factors influencing such a large spread may depend on the cultivation conditions, harvest time, storage method, degree of humidity, etc. [10,11]. The result of the above-mentioned high variability of the properties of fibers is the problem of predicting both the properties of the fibers in various forms and their biocomposites.
In the current reality, modeling the work of a sample made of a specific composite material requires detailed calculations on the micro and macro scale of the material. These problems are not well described in the publications mentioned above. It is important to analyze the properties of the material on a micro scale, choosing an appropriate homogenization model. Then, in the next step, to determine the properties of a representative volume element (RVE), it is possible to simulate the tests for a sample made of a specific composite, taking into account the heterogeneous structure of the material and the filler geometry. Such an advanced type of the analysis may reflect the behavior of the sample during a real test or an operation of a specific product made of a composite reinforced with continuous fibers of a plant origin. A new trend in forecasting and verifying the properties of polymer composites filled with short and continuous fibers of the plant origin may be the use of Digimat software in micromechanical calculations [12,13]. This software uses micromechanical models of Mori-Tanaka and Double Inclusion, based on the Eshelby's solution. Moreover, in the aspect of finite elements, it is possible to generate advanced geometry of reinforcement, reflecting their real structure. An example is the publication [14] where the forecasting of the properties of flax fabric-epoxy resin composite was performed using the Mori-Tanaka and Double Inclusion homogenization model in micromechanical modeling. Moreover, the second significant practical problem that was implemented in the study was the assessment of the impact of the size of the RVE on the received results. The calculations were conducted for the composites with a polymer reinforced with flax fabrics type of: twill 2x2, twill 3x1 and plain. It was noted that the received results depend on: a properly selected homogenization method, the dimension of the RVE, the amount of yarns in the RVE and the sort of weave in the fabric. The received data may be used, among others, to carry out the practical calculations of the work of a real product made of the analyzed composite.
Due to the relatively high variability of the properties of natural fibers, the use of appropriate computational methods to predict and verify the properties of the resulting composites filled with natural fibers constitutes a great problem. This explains a very small number of publications, especially those concerning natural continuous fiber as reinforcement in the polymer matrix. In the work [15], the numerical calculations for simulating the mechanical properties of polymer composites reinforced with long flax fibers in the form of a fabric were presented. By means of a Python script within Salome-Meca 2015 software, the geometry of composite was generated. The bidirectional flax fabric was created and discretized by finite elements type of 4-node tetrahedron. For defects in fibers and interfacial zones of fiber yarns, the law of brittle material was selected. Taking into account the isotropic hardening law and non-local continuum failure mechanics, flax/epoxy composite including fabrics, was modeled with a nonlinear plasticity model. The nonlinear system of equation was solved by using the Newtone-Raphson method. The results of numerical simulations were compared with the experimental data. The results of this paper show that the simulation can properly capture the main failure mechanisms of composites, such as fiber breakage, damage to the polymer matrix and fiber separation at the matrix/fiber connection. Moreover, the numerical method results present a good agreement with the experimental results in terms of elastic and non-linear properties in the stress-strain behavior of polymer composites reinforced with fibrous fabric. The paper [16] analyzed the tensile behavior of epoxy resin composite reinforced with jute fibers using the Siemens PLM NX 10.0 software. This type of materials must have sufficient stiffness to withstand deformation of the geometry under normal loading circumstances. Therefore, the work focused on the analysis of the load -deformation behavior of the prepared jute/epoxy composite. A three-dimensional model of the tested detail, including the arrangement of layers was created. The values of stress-strain characteristics were verified on the basis of the literature. The applied model enabled to achieve good compliance between the simulation and experimental results. It should be noted that the study investigated the influence on the type of fiber, fiber content, layers number, layers angle and fiber shape in the weft on the laminate stiffness. In work [17] the finite element analysis of multi panel structures made of composites reinforced with coir and jute fibers was conducted. The composites were interpreted as isotropic material with uniform Young's modulus and Poisson's ratio. In the Abacus software, the mechanical behavior of a two-panel plate and a six panel box structure under various loading conditions was analyzed. In order to examine the relationships between various parameters, deflections of the panels and stress distributions, exhaustive parametric studies were also performed. The results were compared with the experiment, where a good level of agreement was found. The research [18] presents the simulation of a composite consisted of unidirectional flax in a polypropylene matrix under quasi-static tensile loading, where an anisotropic continuum failure modeling approach was applied to damage a model. An anisotropic strain-dependent material failure model was calculated and applied, based on the experimental damage and strength values. Numerical simulation was performed in the Abacus software, where geometric models of the fabric including yarn geometry were prepared. A comparative analysis of the simulation with the experiment was performed in the field of the stress-strain characteristic, determining the failure strength of the material, where the differences in the values amounted to approx. 3.1%. In paper [19], the mechanical behavior of unidirectional bamboo fibercomposites exanimated in the uniaxial tensile test was predicted. To this end, multi-scale approach calculations were conducted based on the micromechanical analysis of the failure at the matrix and fiber scale. The above-mentioned composite was analyzed under the transverse matrix cracking and the conditions of fiber cracking, taking into account longitudinal splitting resulting from a poor interfacial bonding between the matrix and filler. The improved Weibull model was selected to present the fiber strength distributions that are affected by changing the diameter of the fiber. In order to characterize the progressive damage of the composite, a simulation package was performed using the Monte Carlo method. The accuracy of the selected model to predict the strength of the analyzed composite was noted, based on the comparison of the results from the simulation and the experiment. In work [20], panels consisted of unidirectional flax fibers tapes with polypropylene films were experimentally and simulating studied in terms of the bending and tensile properties. The changing in the properties of the above-mentioned composite was found to be relatively moderate as compared with that of single natural filler. Multi-scale finite element analysis (FEA) calculations for the progressive failure prediction of composite were conducted. The first step was a micromechanical analysis using a representative volume element (RVE) to predict the property of the unidirectional flax layer. Then, using the obtained results, the macro-scale forecasting of the composite layer package properties was carried out. The prepared multi-scale finite element model allowed forecasting the composite bending behavior and the tensile strength, by taking into account the main model determining the damage of the PP-flax fiber composite.
On the basis of the possibilities of micromechanical calculations with the use of homogenization models, it is possible to obtain, among others, the data such as engineering constants and stiffness matrices. These data can be used, inter alia, to define the material data in the programs used to simulate the work of specific utility products and intended for testing. Due to the small number of scientific reports on the possibility of forecasting the properties and performance of products made of composite materials reinforced with plant-based fiber fabrics, simulation studies were carried out to perform micromechanical calculations and then the obtained data were used to simulate the work/load of a specific detail made of the above-mentioned composites. Therefore, it is worth noting that the performed calculations and simulations refer to the actual results of experimental tests in order to verify the obtained results; therefore, the simulation tests performed concerned the properties of composites and products made of them, tested in the paper [21].

MICROMECHANICAL MODELING OF COMPOSITES PROPERTIES
This study was based on the results of experimental work [21] performed on the composites based on natural fabrics (Biotex Flax 400g/cm 2 and Biotex Jute 400 g/cm 2 ) and the Kinetix R240 epoxy resin made by hand lamination.
In order to create a correct numerical model of a composite material based on layers made of a fabric on a resin matrix, a correct yarn model and then a fabric model should be created, followed by a composite model in the last stage. Such a multistage process of building a numerical model allows obtaining a wide range of information both on the geometry and mechanical properties of the material under study. It also allows one to control the properties at each stage (fiber properties, yarn properties, fabric properties and finally the properties of the composite). On the basis of the manufacturer's data, the mechanical properties of the yarns (flax yarn and jute yarn) ( Table 1 and 2) were calculated using the data on the mechanical properties of individual flax and jute fibers ( Table 3).
The next stage of calculations was to obtain the parameters of the fabric with a given weave pattern based on the properties of yarn after homogenization. It was assumed that the distribution of reinforcement in the composite was regular; therefore, it was possible to determine the RVE. The calculations for various RVE sizes using Digimat MF software were performed.
The Double Inclusion first order homogenization scheme was used in the calculations. The RVE unit cell comprised three warp strands and three weft strands. The dimensions of the unit cell were selected on the basis of the research presented in the paper [14]. Figure 1 shows a view of the actual fabrics and the geometric model. At this stage, the data on the properties of the tested fabrics were obtained, which are presented in Table 4.
After the homogenization process, the detailed results were obtained in the form of engineering constants necessary for further calculations regarding compression of composite profiles with a given layer arrangement ( Table 5). The stiffness matrices for the tested composite materials were also obtained ( Table 6 and 7). The received results of the mechanical properties of the composite were compared with those obtained experimentally (Table 8) [21]. Note: E -Young's modulus, f uT -tensile strength, f uC -compressive strength.

NUMERICAL MODELLING
The fastest method to check the loss of stability of a section is the Linear Buckling Analysis. However, the method has limitations. This analysis ignores key factors such as material and geometric nonlinearity problems and imperfections. Therefore, a more demanding but detailed nonlinear analysis was performed in the Ansys Workbench software, with little information about the material tested and its behavior. In order to check the correctness of numerical model with experiment, the geometry and conditions for modelling were based on the tests presented in research [21]. The finite element 3D composite model was designed by using ACP (Pre) module (Fig. 2). In this module, the fabrics and ply lay-out were defined. It was assumed that the thickness of each ply was constant and equaled 0.7 mm. A numerical analysis was carried out for channels with 2, 3, and 4 number of plies for one direction ply angle (0°) and for channels with 2 and 3 number of plies for different ply angle. The finite element model consists of hexahedral elements (HEX20) with maximal edge size of 2 mm. The size of finite elements was defined with convergence analysis under compression conditions with small deformation, where reaction force and maximum stress value was taken into account. The number of elements was limited  Table 6. Stiffness matrix of Biotex Jute 400 g/ m 2 2x2 twill fabric reinforced composite  Table 7. Stiffness matrix of Biotex Flax 400 g/ m 2 2x2 twill fabric reinforced composite  Note: E C -Young's modulus (compression); E T -Young's modulus (tensile); f uC -compressive strength; f uT -tensile strength; V f -volume fraction; ε uT -tensile elongation. by computer calculation power and time. Thus, number of elements is 29400 for 2 plies, 42300 for 3 plies and 55200 for 4 plies respectively.
In numerical model additional plates were considered symmetrically on both ends of the channel (Fig. 2). The contact conditions between the channel and the plate were defined as "bonded". For this type of connection no separation in normal direction and no sliding in tangential direction are possible in contact pairs.
The lower plate was fixed, while for the upper plate displacement in longitudinal direction of 2 mm was established. The material model of composites was assumed as orthotropic elastic type and resin was modelled as isotropic elastic model. The material data were provided by Digimat MF software. For the damage analysis Tsai-Wu 3D criterion was chosen.
An application of orthotropic material model affects the material properties dependence on the direction. Ansys ACP (Pre) allows generating a graphic presentation of material changes in polar coordinates (Fig. 3). It is evident from these diagrams that the properties composite with flax woven fabric shows greater sensitivity on the direction changes than jute woven fabric composite. However, in both cases this variation of properties is small. For elastic modulus, the maximum change is 7.81% for flax, and 3.55% for jute. In the case of shear modulus, the maximum variation is 0.23% for flax and 0.25% for jute, respectively.
The channel test results for one directional plies angle are presented in table 9 and compression curves are presented in Figure 4. The channels compression response is similar to each considered cases. Two stages of compression can be distinguished: the linear response until bulking inflection point is reached and the short nonlinear response to obtain the maximum force for the displacement of 2 mm. In the first stage, the numerical results are close to the experimental data (Fig 4). After bulking initiation, the curves are diverging. No material damaged effect was implemented in the numerical model, thus ultimate force was not reached and material weakening effect was not included. Only failure initiation force was determined by using Tsai-Wu failure criterion ( Table 9). The failure initiation force is obviously lower than ultimate force of channel, thus cannot be used directly for defining ultimate strength of analyzed object. However it allows determining the point at which microcars of material can occur and predict the construction failure.
The numerical results of bulking force are close to the experimental data. The highest differences were observed for 2-layer flax and 3-layer jute, and are 30% and 25%, respectively, in comparison to experimental results. The mean deviation is under 13%. This divergence can be explained by assuming of constant ply thickness in numerical model, while experimental research presents different thicknesses of prepared channels for each case [21]. Two buckle modes were obtained in simulations, which composed of two or three half waves of buckling (Fig. 5). However, no dependency with material or number of layers was found. The same buckling half waves were observed for web and flanges. The results of the numerical analyses are in line with the experimental studies [21].
The results of compression test for channels with different ply angle are presented in Table 10.
The changes of the ply angle for 2-layer channels affect the value of buckling force and failure initiation force (Fig. 6). The maximum changes of buckling force for the flax composite do not exceed 28% and 25% for the jute composite. No significant relationship between ply angle and buckle force and failure initiation force is observed for 3-layer sections (Fig. 7). Maximum difference: in buckling force does not exceed 5%, in failure initiation force it is within range 4-6%. The effect of bulking mode on bulking force and failure initiation force can be found for two layers channels (Fig. 6). For two half waves mode, the buckle force is reduced by approx. 22% for the flax woven fabric, and 16% for the jute woven fabric (Fig 6a), while the failure initiation force is increased by approx. 7% in comparison with three waves buckle mode (Fig 6b). Three layers channels do not show a clear relationship between the buckling mode and the buckling force.

CONCLUSIONS
The procedure for numerical modelling of composites with natural fiber by using Digimat MF and Ansys software was presented. Digimat   MF software allowed determining the mechanical properties of the composite materials with natural woven fabric for numerical modeling based on the properties of the fibers and the matrix. Double inclusion homogenization model was applied caused by reinforcement nature of chosen fabric weave. As a result, it was possible to obtain the information about the mechanical properties of the tested composites in three directions. The calculated values of Young's modulus obtained in the direction of E1 in comparison to the experiment are higher for Biotex Flax by 1.7%, for Biotex Jute by 0.23%. Buckling forces value deviation did not exceed of 13% in comparison to the experiment. The values of force failure initiation are approx. 19% lesser than experimental ultimate forces and only indicate the starting point of material failure. The shape of compression curves received from simulations are close to experiment, until buckling force is reached. Further differences are observed because no material damage was applied in the model. Two types of C -section buckling modes were obtained and they consisted of two or three half-waves of buckling, which is consistent with the experiment. The buckling mode shape effect on the buckle force and failure initiation force for two layers channels is observed. It was noticed that for 2 -layers composites the ply angle effect on buckling force changing in maximum value of 25%. In the composites with 3 layers, the ply angle changes do not have a significant influence (less than 5%) on buckling and failure initiation force.