Abstract
The paper describes an algorithm for numerical estimation of effective mechanical properties in two-dimensional case, considering finite strains. The algorithm is based on consecutive application of different boundary conditions to representative surface elements (RSEs) in terms of displacements, solution of elastic boundary value problem for each case, and averaging the stress field obtained. Effective properties are estimated as a quadratic dependence of the second Piola-Kirchhoff stress tensor upon the Green strain tensor. The results of numerical estimation of effective mechanical properties of plexiglas, reinforced with steel wire, are presented at finite strains. Numerical calculations were performed with the help of CAE Fidesys using the finite element method. The dependence of the effective properties of reinforced plexiglas upon the concentration of wires and the shape of wire cross section is investigated. In particular, it was found that the aspect ratio of reinforcing wire cross section has the most significant impact on effective moduli characterizing the material properties in the direction of larger side of the cross section. The obtained results allow one to estimate the influence of nonlinear effects upon the mechanical properties of the composite.
1. Introduction
The widespread use of plexiglas in engineering applications is due to its properties such as transparency, strength, flexibility, lightness, cheapness, and nontoxicity. In some applications the combination of plexiglas properties of high rigidity and resistance to mechanical and thermal effects is required. In order to increase stiffness and thermal conductivity of the plexiglas and to prevent its spillage due to mechanical or other influences a reinforcement wirework (usually made of steel) is used. Reinforced plexiglas is a composite material. Properties of the material obtained depend primarily upon the material and shape of reinforcing wire. At that, the question that has to be answered is as follows: how to evaluate mechanical properties of reinforced plexiglas while mechanical properties of plexiglas and wires are known provided the shape of wire cross section is known?
The averaging of heterogeneous materials properties has been of interest since the middle of the last century. Theoretical principles of such averaging are described in [1]; in particular, a concept of representative volume element is explained in detail. The studies of that time concerned the effective properties of composite materials in the linear form suitable for description of composites behavior under small strains.
For composites with relatively small volumetric content of filler in the matrix, the existing bilateral assessment of Hashin and Shtrikman [2, 3] is valid; it gives the minimum and maximum values for the compression bulk modulus and a shear modulus of composite material (provided filler concentration and moduli of filler and matrix are known). Mori and Tanaka method [4] is an application of Hashin and Shtrikman conditions to fiber-reinforced material with a continuous matrix in which the volume fraction of the filler particles is low, and the particles shape is mostly spherical. The book of Christensen [5] contains some analytical formulas for evaluating of effective elastic properties of fiber-reinforced, fibrous, and laminated composites and plastic and viscoelastic effects, and effective thermal properties of composites are studied.
At present, the more topical problem is evaluation of effective mechanical properties of heterogeneous materials in nonlinear form, with the help of which it is possible to describe the behavior of composites under finite strains. In the work [6], elastic and plastic properties of the material containing the dispersed microdefects of different orientation in space are studied. In articles [7, 8], effective elastic properties of solids containing cavities of various shapes and orientations in space are estimated. In [9], for the construction of effective constitutive relations of viscoelastic material with periodic structure the finite element method is used; it serves to solve the two-dimensional problem of elasticity theory for a representative volume element, and then the results are averaged. Articles [10, 11] describe how to apply the variational principle to estimate the effective characteristics of multicomponent composites in the form of the strain energy density. Work [12] presents a method of effective characteristics estimation of composite materials in nonlinear form, based on the Tsukrov et al. [7, 8] principle. Work [13] describes a method for constructing nonlinear thermoviscoelastic effective constitutive relations for composites of a periodic structure with a rubber-like matrix. Work [14] describes the application of probability theory methods for estimating of effective mechanical properties of composites having irregular structure. Work [15] presents a method suitable for evaluating both elastic properties and thermal and electrical conductivity of fiber-reinforced materials (the article compares the influence of different parameters of reinforcing particles upon the effective elastic properties and effective thermal and electrical conductivity). In [16], an averaging of elastic properties of inhomogeneous material is performed under nonperiodic boundary conditions considering geometric nonlinearity; and practical implementation is carried out using the finite element method (for two-dimensional case). In [17], this approach is applied to the multiscale case. Article [18] compares different methods for properties averaging of both linear viscoelastic and nonlinear viscoplastic composite materials.
This work proposes an algorithm for construction of effective constitutive relations of composite materials in nonlinear elastic form under finite strains. An evaluation of the effective characteristics is discussed for two-dimensional case, which is a simplification as compared to the three-dimensional case for which the algorithm is described in [19, 20] in detail. Basic principles of this approach are as follows [21–23] (as for two-dimensional case). A rectangular RSE is taken. An effective material is a solid homogeneous material which meets the following condition: if this homogeneous material fills RSE and the source composite fills the exact same RSE, then the average stress over the area in the original composite and the effective homogeneous material are equal under the identical displacements of boundaries of the RSE. Constitutive relations for the effective material (i.e., the effective properties) are plotted as a quadratic dependence of the second Piola-Kirchhoff stress tensor upon the Green strain tensor. To calculate coefficients for this dependence, several sequences of boundary value problems of the nonlinear elasticity should be solved for RSE with given displacements of boundaries.
The paper presents the results of two-dimensional calculations of reinforced plexiglas effective characteristics. Calculations are performed with the help of CAE Fidesys using the finite element method. Plexiglas and its reinforcing steel wire were modeled by the Murnaghan material. A dependence of the effective elastic moduli on the concentration of the reinforcing wire and the shape of its cross section was investigated.
2. Materials and Methods
2.1. Algorithm of Numerical Evaluation of Effective Mechanical Properties of Composite in Two-Dimensional Case at Finite Strains
Let us present the basic designations and relations of the nonlinear theory of elasticity [24] which will be used in describing the algorithm of numerical evaluation of effective mechanical properties of composite material in two-dimensional case at finite strains: , : radius vector of a particle in initial and current states; − : displacement vector; : material coordinates of a particle; : spatial coordinates of a particle; : basis vectors of reference system; = , : basis vectors in initial and current states; = , = : gradient operators in initial and current states; : identity tensor; = = + : deformation gradient; = : the Green strain tensor; : the Almansi strain tensor; : relative change of volume; : sign of transposition; : true stress tensor; : the first Piola-Kirchhoff stress tensor; = : the second Piola-Kirchhoff stress tensor; , : area of RSE in initial and current states; , : boundary of the RSE in initial and current states; , : normal to the boundary in initial and current states.
Let us call the effective (averaged) material such a homogeneous material that will meet the following condition: if this homogeneous material fills RSE and the source composite fills the exact same RSE, then average (over the area) stresses in the original composite and in the effective homogeneous material are equal under identical displacement of faces. Let us call mechanical properties of this material as effective properties.
Using the above definitions and designation, let us describe the algorithm for estimating effective characteristic of composite material for two-dimensional case in the nonlinear form under finite strains.
For RSE in the initial state (before deformation), we will solve a certain number of sequences of boundary value problems of nonlinear theory of elasticity [24]: with boundary conditions
Each sequence of problems corresponds to a certain type of boundary conditions (2) and a certain type of effective strain tensor over the RSE (i.e., to a certain type of effective deformation gradient ). Different problems within the same sequence differ by strain values.
Solving each problem of each sequence we will find the field of the true stress tensor . With the knowledge of it we will calculate the effective stress tensor by averaging of over the area:
The second equality in (3) was obtained with the help of divergence theorem and the fact that the RSE is in equilibrium:
Knowing the deformation gradient that was specified in (2), it is possible to calculate the effective Green strain tensor using the following formula:
In a linear case, the effective constitutive relations are estimated as a linear dependence of true stress tensor upon strain tensor:
In the nonlinear case, the obtained effective true stress tensor is used first for calculation of the effective second Piola-Kirchhoff stress tensor :
In this case, the effective properties are estimated as quadratic dependence of the second Piola-Kirchhoff stress tensor upon the Green strain tensor :
Thus, estimation of the effective properties of the composite in linear case reduces to the calculation of coefficients (6), and estimation of effective properties in nonlinear case reduces to calculation of coefficients and (8).
At that, the following conditions of symmetry are in force for tensor components in (6):(1) (due to symmetry of stress tensor);(2) (due to symmetry of strain tensor);(3) (due to existence of strain energy density).
The same conditions are in force for in (8). For in (8) the following conditions are in force:(1) = (due to symmetry of strain tensor);(2) = (due to symmetry of strain tensor);(3) = (due to symmetry of multiplication);(4) = (due to symmetry of the Piola-Kirchhoff tensor).
Considering the above conditions of symmetry, there are 21 independent constants or and 126 independent constants .
2.2. Finite Element Implementation of the Algorithm
Implementation of the algorithm relies on the use of finite element method [25, 26] for the calculations. For numerical estimation of effective properties, a geometric model of RSE of a composite material is required, which should be “cut out” in the shape of a rectangle whose edges should be parallel to the coordinate planes. The finite element mesh, which is an input for the algorithm, is constructed on this geometric model. Further actions are performed in the Cartesian coordinate system, coordinate planes of which are parallel to the edge of RSE. In this coordinate system, the effective mechanical characteristics will be evaluated.
The algorithm consists of the following logical blocks.
(1) Preparing the Mesh for Calculation. The finite element mesh includes an array of node numbers with their coordinates and an array of numbers of finite elements with the numbers of their constituent nodes. In order to determine the overall dimensions of the RSE, a calculation loop over all nodes of the finite element mesh is conducted, which defines the maximum and minimum abscissas and ordinates of mesh nodes: , , , and . Then the coordinates of the RSE center are calculated:
After that, the model is shifted to put its center to the origin (i.e., the coordinates of all mesh nodes are decreased by and ). Also, overall dimensions of the model are determined by formulas:
In addition, the accuracy for determination of boundary nodes and edges is calculated for the same block. For this purpose a calculation loop over all finite elements of the mesh is conducted, which defines the minimum area of an element . The accuracy is calculated by formula
The value of naturally depends on how fine mesh is used for calculation.
(2) Formation of Lists of Boundary Nodes. For further calculations it is necessary to form a list of nodes on boundaries of the RSE. Particular lists depend on type of boundary conditions applied.
If the problem is solved for nonperiodic boundary conditions, all nodes located on the boundaries are put to common list. For this purpose the calculation loop is conducted over all mesh nodes, at each step of which the following condition is checked: OR OR OR . Here, and are coordinates of the node. If the condition is true, this means that the node is on the boundary of the RSE, and its number is added to the list of boundary nodes.
If the problem is solved for periodic boundary conditions, four lists of numbers of nodes are created (according to number of edges). Similarly, the calculation loop is conducted over all mesh nodes, at each step of which the following condition is checked:(1);(2);(3);(4).
If the condition with number is true for a node, the node number is added to the list with the number .
Then, a pair of corresponding nodes (i.e., those, the projection of which is closest to each other) are formed from opposite edges of the RSE. The pairs are recorded to two lists.
For this purpose, a pass through all nodes of the edge is conducted. For each such node a node is found: the closest one to the point on the edge that is opposite to the node . In other words, for the node with coordinates () (which lies on the edge ) the closest one (from nodes lying on the edge ) to the point with coordinates () is searched for. When such a node is found, the pair () is added to the first list of pairs of nodes. Then the loop is conducted over all nodes of the edge : if the node of this edge is already the second element of one of the pairs, we can proceed to the next one; but if not, then you need to find its corresponding node on the edge and add it to the first list of pairs.
Similarly, a pass through all nodes of the edge is conducted. For each such node a node is found: the closest one to the point on the edge that is opposite to the node . In other words, for the node with coordinates () (which lies on the edge ) the closest one (from nodes lying on the edge ) to the point with coordinates () is searched for. When such a node is found, the pair () is added to the second list of pairs of nodes. Then the loop is conducted over all nodes of the edge : if the node of this edge is already the second element of one of the pairs, we can proceed to the next one; but if not, then you need to find its corresponding node on the edge and add it to the second list of pairs.
(3) Formation of Lists of Boundary Edges. For further averaging over an edge it is necessary to form a list of edges of finite elements on the boundaries of the RSE. Since the finite element mesh data do not usually contain an array of edges of elements, four lists of pairs are formed. For this purpose a calculation loop over all finite elements of the mesh is conducted, and a loop over all edges of the element inside it is conducted, too. If all nodes of (local) element edges lie on one of the edges of the RSE (global), a pair is added to an appropriate list. This takes into account that the elements are of different geometric shape and different order of approximation; that is, different elements have different number of edges and different number of nodes on them.
(4) Application of Boundary Conditions. For RSE we will solve a certain six sequences of boundary value problems of nonlinear theory of elasticity [24]. Different problem sequences differ by types of applied boundary conditions (i.e., type of effective strain tensor over the RSE). Different problems within the same sequence differ in the amount of strain while being of the same type. The following types of strains are applied:(1): strain or compression along the axis ;(2): strain or compression along the axis ;(3): shear in plane;(4): composition of tensions or compressions along two axes: and ;(5); : composition of tensions or compressions along the axis and shear in plane;(6); : composition of tensions or compressions along the axis and shear in plane,
where is the amount of strain.
There are 6 sequences, each of which contains 3-4 problems or more, if required. Different problems within the same sequence differ in the amount of strain .
On the basis of effective strain tensor, the effective deformation gradient is calculated by formula (5). Since the strain tensor is symmetric and deformation gradient in common case is asymmetric, for certainty it is assumed that the gradient is upper triangular. Then (5) can be figured out by components in the following form:
Formulas for gradient’s components in an explicit form are as follows:
For nonperiodic boundary conditions there is one common list of numbers of all nodes on the boundary of the RSE. Boundary condition (2) is applied to each node in this list in the form of rigidly given displacements in all two directions.
For periodic boundary conditions there are two lists of pairs of nodes on the boundary of the RSE. For each pair of nodes from the first list () the relationship with displacement is specified in the following form:
For each pair of nodes from the second list () the relationship with displacement is specified in the following form:
(5) Solution of Boundary Value Problem of Elasticity Theory and Averaging of Results. After application of boundary conditions to the RSE, the actual numerical solution of each boundary value problem of the elasticity theory of each sequence is performed using the finite element method. When the solution is found, it is necessary to average obtained field of stress tensor over the area.
First, the area of the model in the final state (i.e., after deformation) is calculated using the formula:
In other words, the area of model in final state represents the area in initial state multiplied by the determinant of effective deformation gradient.
Stress tensor is averaged over the area using formula (3). Since the integration is performed on the finite element model, a transition from integral over the entire boundary to the sum of the integrals is done over all boundary edges:
Hereby one can calculate the effective true stresses tensor in the RSE. Knowing it, it is possible to calculate the effective second Piola-Kirchhoff stress tensor using formula (7).
Blocks 4 and 5 of the algorithm should be run in double loop: by type of strain (1 to 6) and by amount of strain (from 1 to 3…4 or more). The result is as follows: for each boundary value problem for each sequence an effective Green strain tensor was set, and as a result the effective second Piola-Kirchhoff stress tensor was obtained. The resulting Piola-Kirchhoff tensors are stored for each problem for later calculation of effective properties in the form of their relations.
(6) Plotting of Piola-Kirchhoff Tensor versus Green Tensor Using the Least Squares Method. Since the problems within one sequence differ by the amount of strain only, for each problem of each sequence the dependence of the effective second stress Piola-Kirchhoff tensor versus characteristic amount of strain is plotted:
The dependence is constructed using the least squares method, which allows calculating coefficients and .
(7) Calculation of Effective Elastic Moduli of First and Second Orders. In block 6 for each problem of each sequence, we have constructed the dependence (18) of the effective second Piola-Kirchhoff stress tensor versus the characteristic amount of strain, . Effective properties are evaluated as (8).
The relationship between the calculated coefficients and in (18) for each problem sequence and desired effective elastic moduli (of the first order) and (of the second order) is as follows:(1) = ⇒ = + = + ;(2) = ⇒ = + = + ;(3) = ⇒ = + = + ;(4) = ⇒ = = + ;(5) = ⇒ = + = + ;(6) = ⇒ = + = + .
Formulas for calculation of coefficients in an explicit form are as follows:(1) = ;(2) = ;(3) = .
For calculation of coefficients it is necessary to solve the system of 6 linear algebraic equations:(1) = ;(2) = ;(3) = ;(4) = ;(5) = ;(6) = .
This system is solved analytically, and the solution is as follows: = ; = ; = ; = ; = ; =
Thus, the described algorithm allows performing numerical (using the finite element analysis) estimation of effective elastic moduli (of the first and the second order) for the composite material in the two-dimensional case. On the basis of the described algorithm we have developed a software module Fidesys Composite as part of CAE system Fidesys [27] intended for numerical estimation of effective mechanical properties of composite materials.
The proposed algorithm assumes that the computations for each geometrical pattern are performed separately. For composites of periodical structure it is sufficient to perform computations once. For composites of random structure with statistically uniform distribution of inclusions one can use the approach that is based on the ensemble averaging over a fixed number of configurations [21, 22]. The additional averaging over all possible directions in the plane of deformation [21, 22] can be performed analytically and results in effective constitutive equations for transversely isotropic materials. These methods of averaging permit one to avoid the huge amount of computations for irregular composites.
3. Results
With the help of the developed module Fidesys Composite, two series of finite element analyses of effective characteristics of plexiglas reinforced with steel wire were conducted for two-dimensional case in the nonlinear form under finite strains. A rectangular RSE of plexiglas with a rectangle of steel modeling wire was considered in all calculations. A level of strains applied to the RSE was 1%, 2%, and 3%. Mechanical properties of plexiglas and steel were described using the Murnaghan constitutive relations [24]:
The constants = 1.09 · 105 MPa, = 0.818 · 105 MPa, = −0.29 · 105 MPa, = −2,4 · 105 MPa, and = −2.25 · 105 MPa [24] were used for steel, and the constants = 0,39 · 105 MPa, = 0.186 · 105 MPa, = −0.013 · 105 MPa, = −0.07 · 105 MPa, and = 0.063 · 105 MPa [24] were used for plexiglas.
3.1. Dependence of Effective Properties versus Concentration of Wires
Dependence of effective properties of reinforced plexiglas versus the concentration of the reinforcing wires was studied. RSE was a rectangle 10 × 5 mm. Cross section of wire is square; size of square side ranged from 0.25 mm to 3 mm.
Let us see graphs for linear coefficients , , and (see Figures 1–4).
The graphs show that the coefficient defining the material behavior during the deformation along the axis almost linearly depends on concentration of wires (Figure 1): the greater the concentration, the more . The dependence of coefficient responsible for the behavior of the material during the tension along the axis versus the concentration is similar (Figure 4). Similarly, coefficients (Figure 3) defining the behavior of the material under shear strain and coefficient (Figure 2) are dependent upon concentration of wires.
Let us see graphs for nonlinear coefficients , , and as well (see Figures 5–7).
As can be seen from the graphs, the coefficient (which determines the nonlinearity of the material during tension along the axis ) increases in modulus with increasing of concentration of wires (Figure 5). The coefficient (responsible for the nonlinearity during tension along the axis ) increases in modulus even stronger, and the dependence of this coefficient on concentration of wires is clearly nonlinear (Figure 7): this is probably due to the fact that the width of the RSE is half as much as length. Other nonlinear coefficients (e.g., , the graph for which is shown in Figure 6) depend on concentration significantly lesser.
3.2. Dependence of Effective Properties versus the Form of Cross Section
Dependence of effective properties of reinforced plexiglas versus the form of cross section of the reinforcing wires was studied. RSE was a square 10 × 10 mm. The rectangular form of wire cross section was considered, the ratio of rectangle sides varied from 1 : 1 to 1 : 10. At that, the cross-section area of the wire was 9 mm2.
Let us see graphs for linear coefficients and (see Figures 8 and 9).
The graphs show that the coefficient grows monotonically and almost linearly as the cross section of the reinforcing wire is “pulled” along the axis (Figure 8). At that, the coefficient decreases (Figure 9) monotonically but nonlinearly: at first, this coefficient falls sharply (with the change of aspect ratio from 1 : 1 to 1 : 3) and then changes slightly. Other linear coefficients are slightly dependent upon the shape of wires section as shown by calculations.
Let us see graphs for nonlinear coefficients , , and (see Figures 10–12).
As can be seen from the graphs, the coefficient increases in modulus monotonically as you “pull” the cross section of the wires along the axis , (Figure 10), and the dependence is close to linear one. At that, the module of coefficient decreases (Figure 12) also monotonically but nonlinearly: similar to the coefficient , the coefficient decreases in modulus strongly with decreasing of aspect ratio of the rectangle from 1 : 1 to 1 : 3; then it changes weakly. The graph for the coefficient , defining the nonlinearity of material under shear strains, looks more complicated (Figure 11): first, upon the increase of the aspect ratio of the rectangular cross section, the modulus of this coefficient drops by approx. six times, and then it rises again by seven times. The minimum of modulus is observed when the aspect ratio is approximately 1 : 6.
The points of inflection in the figures showing stiffness constants as a function of the width/length ratio of the inclusions arise apparently due to computational errors.
3.3. Comparison with Results Obtained by Other Methods
Within the framework of linear elasticity and small strains, the obtained results are compared with the results of Christensen [5]. The materials of matrix and fibers (wires) were assumed to be isotropic, and the comparison was performed for the plane strain. In our computations RSE was a square, and the circular inclusion was in the center of RSE. The concentration of the reinforcing wire was 10%. The mechanical properties of inclusions (wires) are Young’s modulus = 2000 MPa and Poisson’s ratio = 0.2. The mechanical properties of a matrix are Young’s modulus = 2 MPa and Poisson’s ratio = 0.3. The effective elastic modulus is given in Table 1.
One can see from Table 1 that the difference between the numerical solution and the analytical results given in [5] does not exceed 3%.
Another comparison was performed for the case of finite strains. The pure shear of a fiber-reinforced elastomeric composite with rigid circular fibers (wires) was considered. The matrix material was incompressible, and the mechanical properties of the matrix material were described by the neo-Hookean potential. The stresses were averaged over the RSE, and the first Piola shear stress was computed. The results were compared with the estimations of Avazmohammadi and Ponte Castañeda [28]. These results are shown in Table 2. In this table is a shear strain and is the concentration of wires. The results of the comparison may be estimated as satisfactory.
4. Discussion
Thus, the paper presents the algorithm for numerical evaluation of the effective mechanical properties of nonlinear elastic solids in two-dimensional case (under plane strain). The novelty of this algorithm is determined by considering of nonlinear effects. Both physical and geometric nonlinearity is taken into account. Effective constitutive relations are presented in the form of quadratic dependence of averaged strains versus stresses. Determination of effective modules is reduced to solution of sequences of nonlinear boundary value problems of elasticity for loads of different types and sizes. For solving boundary value problems the finite element method is used, it is implemented in CAE system Fidesys. Results of calculation of reinforced plexiglas effective characteristics, presented in the paper, confirm the efficiency of the algorithm. The influence of concentration of reinforcing wires and wire shape upon the effective properties is studied. In particular, it was found that the aspect ratio of reinforcing wire cross section has the most significant impact on effective moduli characterizing the material properties in the direction of larger side of the cross section.
The obtained results allow estimating the influence of nonlinear effects upon the mechanical properties of the composite for different amounts of strain. For example, in Figures 8 and 10 one can see that for elongation of 10% in the direction of the axis the amendment for adoption of nonlinear effects for the stress in the direction of the same axis is 1% for rectangular wires, and for wires cross section with an aspect ratio of 10 : 1 this amendment will amount to ca. 16% under the same strain.
In the future it is planned to perform similar calculations for other composite materials.
Competing Interests
The authors declare no competing interests.
Acknowledgments
The research for this paper was performed in Fidesys LLC and was financially supported by Russian Ministry of Education and Science (Project no. 14.579.21.0076; Project ID RFMEFI57914X0076).