Design of Reinforcement in Nano- and Microcomposites

The application of numerical homogenization and optimization in the design of micro- and nanocomposite reinforcement is presented. The influence of boundary conditions, form of a representative volume element, shape and distribution of reinforcement are distinguished as having the crucial influence on a design of the reinforcement. The paper also shows that, in the optimization problems, the distributions of any design variables can be expressed by n-dimensional curves. It applies not only to the tasks of optimizing the shape of the edge of the structure or its mid-surface but also dimensional optimization or topology/material optimization. It is shown that the design of reinforcement may be conducted in different ways and 2D approaches may be expanding to 3D cases.


Introduction
The group of non-homogeneous materials includes, e.g., composite materials, liquid mixtures, media with liquid and solid particles, as well as porous media. Reflecting deeper on the structure of a material, one can conclude that all materials are heterogeneous on a given scale. However, for some media, inhomogeneity can be observed only after moving to the atomic scale. Heterogeneity affects the physical processes occurring in the selected medium and causes difficulties in the description of such a medium. Due to the above, the idea of determining, if possible, a macroscopically equivalent medium has emerged.
In place of a non-homogeneous medium in nano-or microscale, a homogeneous medium in macro-scale is introduced. The terms of the microscopic and macroscopic medium concern the distinction of the scale of heterogeneity from the scale of homogeneity. The method of theoretically moving from a heterogeneous scale to a homogeneous scale is called homogenization. Generally, the homogenization techniques deliver the properties of composite considering the properties of constitutes and their volume fractions. The homogeneous body with effective material properties is introduced instead of the heterogeneous body [1][2][3]. The two most straightforward and typical models of micromechanics are the parallel and series models yielding the upper and lower bounds of the composite properties. The parallel model introduced by Voigt [4] for the estimation of the average constants for polycrystals is based on a uniform strain and gives the upper-bounds (the known rule of mixture). The assumption of Reuss [5] also known as the series model is based on constant stress and gives the lower-bounds [2]. These approximations consider only the elastic properties of constitutes and their volume fractions. In the homogenization procedure, we can also distinguish direct and indirect approaches. In the direct homogenization, the volume averaging of stress, strain and energy is applied and the effective properties are calculated according to the definition of effective properties of the composite [6]. This approach is applied in numerical (finite element method-FEM) homogenization for calculation of local field quantities, whereas the geometry and properties of microstructure can be arbitrary. In the indirect homogenization, the Eshelby's equivalent inclusion method is used [7]. The effective properties of the composite are computed based on volume fraction and geometry of the reinforcement as well as the properties of components. Hill [8] developed the self-consistent method for a single ellipsoidal inclusion embedded in an infinite medium. The author presented that the stiffness tensor components for composite lie somewhere in the interval between the Reuss and Voigt bonds, regardless of the geometry. Later, Christensen [9] presented the generalized self-consistent scheme. Composite constituents act together, and the material properties are the result of their interaction. Thus, Hashin and Shtrikman [10] proposed how to determine narrower bonds for the homogenized elastic properties for composite with isotropic constituents taking into account also the interaction between the elastic properties. They applied the principle of minimum potential energy and the idea of polarization. Halpin [11] presented Hill's moduli in the more straightforward analytical form known as Halpin-Tsai equations and extended its application to a wide variety of reinforcement geometries [12]. Mori and Tanaka [13] proposed that the average strain in the interacting inclusions can be approximated by that of a single inclusion (fiber) embedded in an infinite matrix subjected to the uniform average matrix strain. This method was next reformulated by Benveniste [14]. In Mori-Tanaka theory, concentration matrices are given by the solution of a single inclusion in an infinite matrix. Taking into account a two-phase composite with an ellipsoidal inclusion, the Eshelby's tensor is applied to the determination of the elastic field in a particle. Explicit expressions of the components of Eshelby's tensor are presented in work by Mura [15]. Both Bensoussan et al. [16] and Suquet [17] discussed an asymptotic homogenization theory being a mathematical description of periodic composites. Aboudi [18] has presented a unified micromechanical theory using interacting periodic cells to predict the overall behavior of composites. The homogenization is also applied in the characterization of nanocomposites [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34].
Using the homogenization theory, we can determine the effective mechanical properties of a representative volume element (RVE) made of fiber and matrix. The shape of the edge of the area occupied by the fiber has a significant effect on the homogenized mechanical properties for the elementary cell. Consequently, it also affects the stiffness of the laminate. Therefore, the design of the composite involves not only the tasks of optimizing the shape of the edge of the structure or its mid-surface but also dimensional optimization or topology/material optimization-see, e.g., references [35][36][37][38]. The functionally gradient material (FGM) distributes the material functions throughout the material body to achieve the maximum mechanical properties [39][40][41]. FGMs are mostly characterized as structures having different properties in the thickness direction; however, other propositions consider density grading in a width or length direction. Moreover, diagonal grading (in two directions) is also discussed.
The application of the optimization in homogenization problems has been introduced by Bendsoe and Kikuchi [42] and later discussed in papers, e.g., Refs. [43][44][45][46][47]. The homogenization method provided much valuable information used not only in the theoretical description but also in engineering practice. The use of homogenization methods helps in the design of micro-and nanocomposite materials with the desired properties.
The use of numerical homogenization and optimization in the modeling of micro-and nanocomposite reinforcement is presented here. From the mathematical and also numerical point of view, the description of microcomposites and nanocomposites is almost identical. The difference is only related to the real dimensions of the analyzed unit cells. The homogenization problem is solved based on the finite element method and compared with other approaches. It is shown that the design of reinforcement may be conducted in different ways and 2D approaches may be expanding to 3D cases.
In this paper, the following issues influencing the reinforcement design in homogenization problem are discussed: 1.
Boundary conditions of the representative volume element, 2.
Form of the representative volume element,

Preliminary Remarks
In homogenization, a heterogeneous body is replaced with a homogeneous one based on physical relationships between stresses and strains- Figure 1. Here, through a heterogeneous body, we understand the system in which there are two different phases, e.g., matrix and reinforcement. In contrast, a homogeneous system is understood as the system that characterizes effective material properties. The effective Hook's law for a composite is defined as: where the overbar denotes the average over the unit cell and and ̅ is the average stress and strain tensor, respectively, [ * ] is the effective elastic modulus, for which a total number of independent components is controlled by the prescribed symmetry, i, j = 1,2,3.
In the micromechanical analysis, the geometry of the arrangement of fibers in the matrix is an essential condition to develop a material model for the analysis. In the present study, it is assumed that the fibers and the matrix are only two phases in the composite, the fibers and matrix are perfectly bonded, fibers are aligned and have a uniform diameter along its length. Therefore, the transverse isotropy is assumed. Both constituents in the composite comply with Hooke's law.
The transversely isotropic material is characterized by a set of five equations having five effective independent stiffness moduli * , * , * , * , * , which can be expressed in terms of five independent engineering constants, such as the axial and transverse Young's moduli * and * , the axial and transverse Poisson's ratios * and * , and the axial shear modulus * -see Ref.
[2]. In the numerical homogenization, the average stresses and strains are calculated using the volume averaging procedure: where N is a total number of integration points in the unit cell, , are stress and strain component, respectively, defined at integration points k and is the connected volume. The effective Hook's law for a composite is defined as:

Boundary Conditions
where the overbar denotes the average over the unit cell and σ ij and ε ij is the average stress and strain tensor, respectively, [C * ] is the effective elastic modulus, for which a total number of independent components is controlled by the prescribed symmetry, i, j = 1,2,3.
In the micromechanical analysis, the geometry of the arrangement of fibers in the matrix is an essential condition to develop a material model for the analysis. In the present study, it is assumed that the fibers and the matrix are only two phases in the composite, the fibers and matrix are perfectly bonded, fibers are aligned and have a uniform diameter along its length. Therefore, the transverse isotropy is assumed. Both constituents in the composite comply with Hooke's law.
The transversely isotropic material is characterized by a set of five equations having five effective independent stiffness moduli C * 11 , C * 12 , C * 22 , C * 23 , C * 44 , which can be expressed in terms of five independent engineering constants, such as the axial and transverse Young's moduli E * A and E * T , the axial and transverse Poisson's ratios ν * A and ν * T , and the axial shear modulus G * A -see Ref. [2]. In the numerical homogenization, the average stresses and strains are calculated using the volume averaging procedure: where N is a total number of integration points in the unit cell, σ k ij , ε k ij are stress and strain component, respectively, defined at integration points k and V k is the connected volume.

Boundary Conditions
The three classical boundary conditions applied in the analysis of RVE include: • Linear displacement boundary condition (Dirichlet condition), • Constant traction boundary condition (Neumann condition), • Periodic boundary condition.
Historically, the micromechanical models are the mechanical or engineering models. In the 1970s, a mathematical equivalent to micromechanical methods appeared, being the asymptotic homogenization theory. The fundamentals of asymptotic theory can be found, e.g., in Sanchez-Palencia [48], Benssousan et al. [16], and Suquet [17], among others. The asymptotic homogenization theory used periodic boundary conditions in the modeling of composites. Suquet [17], among others, showed that the "plane-remains-plane" boundary conditions are over-constrained conditions. It was shown that characteristic modes of deformation do not result in plane boundaries after deformation. Firstly, used "plane-remains-plane" boundary conditions on RVEs are only valid for the symmetric RVE loaded by normal traction, e.g., Xia et al. [49]. Homogenization theory, which implies the periodic boundary conditions, conducts to the more accurate results than the use of homogeneous displacement and homogeneous traction boundary conditions, e.g., Nguyen et al. [50]. Periodic boundary conditions are simplified to ordinary boundary conditions for the symmetric RVEs [6].
Considering the numerical homogenization, the effective strain and tangent operator associated with the strain tensor at any point in the macro scale are computed by the formulation of the boundary problem related to the microscale at that point [50]-see Equations (1)-(3). In the numerical homogenization (finite element method), the constitutive macroscopic assumption is not required, but, in this case, the issue of the appropriate use of periodic boundary conditions emerges, e.g., Sun and Vaidya [51]. The periodic boundary conditions are formulated as constrained equations of displacement differences for each pair of the opposite parallel boundary surfaces of the RVE. The periodic boundary conditions are the most efficient in terms of convergence rate from the three kinds of boundary condition for the growing size of the RVE [50].
In the present studies, we have focused on the numerical homogenization based on the finite element analysis of RVE in 2D and 3D problems. After solving the assumed numerical problem for specific displacement/force boundary conditions, the distribution of stresses and strains at each point of the inhomogeneous two-phase medium are obtained. The stress/strain values are variable, i.e., σ ij = σ ij (x, y, z) and ε ij = ε ij (x, y, z). The homogenization of such a body requires the assumption of a specific material model of an equivalent homogeneous body to determine the required number of equivalent material constants. There is a full analogy here with the methodology of experimental studies of anisotropic bodies. To demonstrate the finite element (FE) approach, we discuss a problem of transversaly isotropic composites-see also Chwał and Muc [31].
The displacements and tractions used at the boundaries of the RVE are generalized by [52]: where S is the bounding surface, x j is the surface coordinate, and n j is the component of the normal vector to S. To obtain the material constants for assumed material model, the appropriate, e.g., displacement u i should be applied on the RVE that would produce uniform stresses σ ij . In numerical analyses, a series of FE simulations were carried out, in which RVEs were subjected to the various boundary displacements. The boundary conditions applied to the faces of the unit cells are presented in Table 1. The characteristic faces of the quarter of RVE are: x 1 = 0, x 1 = a, x 2 = 0, x 2 = b, x 3 = 0 and x 3 = c-see Figure 2. The u 1 , u 2 and u 3 denote the displacements along x 1 , x 2 and x 3 directions, respectively, and e is the assumed small displacement in the elastic regime. The effective material properties were computed using the appropriate boundary conditions (Table 1) and Equations (1)-(3).
Materials 2019, 04 , x FOR PEER REVIEW 5 of 24  Assuming an 2D orthotropic body, there are four independent material constants, i.e., axial and transverse Young's moduli, in-plane Poisson ratio and in-plane shear modulus which are determined experimentally by three types of tests: (i) elongation in a direction parallel to the fiber, (ii) elongation perpendicular to the fiber; and (iii) shear in-plane. In numerical homogenization of such a body, it is necessary to build three independent numerical models that perform elongation in a direction parallel to the fiber, elongation in the direction perpendicular to the fiber and in-plane shear. We determine the average values of stresses and strains in both directions-Equations (2) and (3), and then, using the physical equation, we calculate the necessary material constants.

Form of the Representative Volume Element
The distribution of reinforcement in micro-and nanocomposites is of high attention and has a significant influence on their mechanical behavior. In fibrous composite materials, the fiber distribution is mostly random and usually unknown. Thus, some idealizations of fiber distribution in the polymeric matrix are assumed in the analysis of composite materials.
According to Hill [53], in the perfect distributions of fibers, it is seen that, due to symmetry and periodicity of these arrays, one representative array can be selected to analyze the lamina at the microscale. Thus, from the whole specimen, a representative subregion called representative volume element (RVE) can be selected- Figure 3. Additionally, this RVE as a volume of material statistically characterizes a homogeneous material. The above approach is used in an analysis of microcomposites and nanocomposites.
There are many ways to idealize the cross-section of a lamina [54,55]. The usually preferred arrangements of unidirectional fibers in the matrix are square and hexagonal arrays. Two popular Assuming an 2D orthotropic body, there are four independent material constants, i.e., axial and transverse Young's moduli, in-plane Poisson ratio and in-plane shear modulus which are determined experimentally by three types of tests: (i) elongation in a direction parallel to the fiber, (ii) elongation perpendicular to the fiber; and (iii) shear in-plane. In numerical homogenization of such a body, it is necessary to build three independent numerical models that perform elongation in a direction parallel to the fiber, elongation in the direction perpendicular to the fiber and in-plane shear. We determine the average values of stresses and strains in both directions-Equations (2) and (3), and then, using the physical equation, we calculate the necessary material constants.

Form of the Representative Volume Element
The distribution of reinforcement in micro-and nanocomposites is of high attention and has a significant influence on their mechanical behavior. In fibrous composite materials, the fiber distribution is mostly random and usually unknown. Thus, some idealizations of fiber distribution in the polymeric matrix are assumed in the analysis of composite materials.
According to Hill [53], in the perfect distributions of fibers, it is seen that, due to symmetry and periodicity of these arrays, one representative array can be selected to analyze the lamina at the microscale. Thus, from the whole specimen, a representative subregion called representative volume element (RVE) can be selected- Figure 3. Additionally, this RVE as a volume of material statistically characterizes a homogeneous material. The above approach is used in an analysis of microcomposites and nanocomposites.
There are many ways to idealize the cross-section of a lamina [54,55]. The usually preferred arrangements of unidirectional fibers in the matrix are square and hexagonal arrays. Two popular idealizations are presented in Figure 3a,b. The nine independent constants are presented from the 3D constitutive equations considering unidirectional lamina being the orthotropic material in idealizations are presented in Figure 3a,b. The nine independent constants are presented from the 3D constitutive equations considering unidirectional lamina being the orthotropic material in nature. Additionally, for a transverse isotropic body, there are five independent constants. Here, it is assumed that the isotropic plane (2-3) is perpendicular to the axial direction (1) of the fibrous reinforcement- Figure 3a.
(a) (b) Here, the attention is focused on the characterization of elastic properties of epoxy nanocomposites with a small amount of carbon nanotubes involving FE modeling contrasted with classical micromechanical approaches. The present calculations consider our earlier studies on nanostructures [23,32,[56][57][58][59][60] and nanocomposites [20,24,25,31]. The topic of polymeric nanocomposites with carbon nanotubes (CNTs) is widely discussed. Mostly the volume fraction of CNTs in nanocomposites is low. At higher CNT fractions, the mechanical properties of the nanocomposite were observed to deteriorate due to the formation of CNT agglomerates, which act as stress concentrators. Since the purpose of the analysis is to evaluate the material constants for nanocomposites, the prepared model should be able to simulate the stress-strain behavior. For this work, it was assumed that the reinforcement is in the form of straight aligned carbon nanotubes uniformly distributed in the epoxy matrix- Figure 3. The behavior of the reinforcement and the matrix was assumed to be isotropic, homogenous and linearly elastic. The typical values of stiffness moduli for carbon nanotube and epoxy resin have been used and are listed in Table 2. In the literature, the effect of carbon nanotubes alignment and waviness on the effective properties of nanocomposite is also discussed, see, e.g., the work by Savvas et al. [34]. Table 2. Material properties of carbon nanotube and epoxy resin used in finite element analysis.

Material Constant
Carbon Nanotube Matrix Young's modulus (GPa) 1000 3.2 Poisson's ratio 0.3 0.3 Here, the attention is focused on the characterization of elastic properties of epoxy nanocomposites with a small amount of carbon nanotubes involving FE modeling contrasted with classical micromechanical approaches. The present calculations consider our earlier studies on nanostructures [23,32,[56][57][58][59][60] and nanocomposites [20,24,25,31]. The topic of polymeric nanocomposites with carbon nanotubes (CNTs) is widely discussed. Mostly the volume fraction of CNTs in nanocomposites is low. At higher CNT fractions, the mechanical properties of the nanocomposite were observed to deteriorate due to the formation of CNT agglomerates, which act as stress concentrators. Since the purpose of the analysis is to evaluate the material constants for nanocomposites, the prepared model should be able to simulate the stress-strain behavior. For this work, it was assumed that the reinforcement is in the form of straight aligned carbon nanotubes uniformly distributed in the epoxy matrix- Figure 3. The behavior of the reinforcement and the matrix was assumed to be isotropic, homogenous and linearly elastic. The typical values of stiffness moduli for carbon nanotube and epoxy resin have been used and are listed in Table 2. In the literature, the effect of carbon nanotubes alignment and waviness on the effective properties of nanocomposite is also discussed, see, e.g., the work by Savvas et al. [34]. Table 2. Material properties of carbon nanotube and epoxy resin used in finite element analysis.

Material Constant Carbon Nanotube Matrix
Young's modulus (GPa) To compute the effective properties of CNT/epoxy nanocomposite, 3D quarter unit cells for square and hexagonal arrays were built and analyzed in the ABAQUS package. The typical dimensions for the quarter of RVE with a square array of CNTs are presented in Figure 3. The same characteristic notations have been used for RVE with a hexagonal array, namely: a-1 2 of the length, b-1 2 of the width, c-1 2 of the high of RVE. However, in the numerical analysis of nanocomposite, the nanotube has the circular cross-section with R i -the inner diameter of CNT, and t-the thickness of CNT-see Table 3. Other characteristic parameters of the conducted FEA are presented in Table 3.  To compute the effective properties of CNT/epoxy nanocomposite, 3D quarter unit cells for square and hexagonal arrays were built and analyzed in the ABAQUS package. The typical dimensions for the quarter of RVE with a square array of CNTs are presented in Figure 3. The same characteristic notations have been used for RVE with a hexagonal array, namely: a-½ of the length, b-½ of the width, c-½ of the high of RVE. However, in the numerical analysis of nanocomposite, the nanotube has the circular cross-section with Ri-the inner diameter of CNT, and t-the thickness of CNT-see Table 3. Other characteristic parameters of the conducted FEA are presented in Table 3. Element type C3D8R-8-node linear hexahedrons with reduced integration FE model of RVEs was built using 8-node linear hexahedrons with reduced integration-C3D8R solid finite elements available in the ABAQUS package. A double layer of elements on the CNT was used for both square and hexagonal RVE- Figure 4. The general structural analysis was conducted.  The results of the current calculation of the effective engineering constants of CNT/polymer nanocomposites are listed in Table 4.  The results of the current calculation of the effective engineering constants of CNT/polymer nanocomposites are listed in Table 4. Table 4. Comparison of the effective engineering constants of CNT/polymer nanocomposites. The present results based on the numerical homogenization for the square and hexagonal arrays are contrasted with results from micromechanical models such as Vanin model [61] (see Equation (A1)). The analytical micromechanical computations were conducted based on the boundary conditions presented in Table 1.

For axial Young's modulus E *
A and Poisson's ratio ν * A , FEA results are almost identical for the square and hexagonal arrays and are the same as from the rule of mixture applied in the Vanin model. Other elastic moduli, the transverse Young's modulus E * T , the axial shear modulus G * A and the transverse Poisson's ratio ν * T represent the main challenge for the researchers. Here, the predicted values of E * T , from FE studies are higher than the micromechanical predictions. The current FE model leads also to the relatively high value of the axial shear modulus G * A that is not predicted by micromechanical models. Lower discrepancy is observed for the transverse shear modulus G * T . The transverse Poisson's ratio ν * T calculated numerically for the square and hexagonal arrays show higher values than computed according to micromechanical model.

Optimization Problems with the Use of 2D Curves
In optimization problems, distributions of arbitrary design variables can be expressed through n-dimensional curves. This applies not only to the problems of optimizing the shape of the edge of the structure or its central surface but also dimensional optimization (e.g., thickness) or topology/material optimization (distributions: density, the volume fraction of strengthening in the material, the material with different Young's modulus)-see Figure 5. In the latter case, the curves form the boundary between materials with different physical and chemical properties. In cases of simultaneous analysis of optimization, e.g., shape and dimensional, we introduce two types of curves: one describes the variations of the shape of the middle surface, and the other changes the thickness-see Figure 5e. The present results based on the numerical homogenization for the square and hexagonal arrays are contrasted with results from micromechanical models such as Vanin model [61] (see Equation (A1)). The analytical micromechanical computations were conducted based on the boundary conditions presented in Table 1.
For axial Young's modulus * and Poisson's ratio * , FEA results are almost identical for the square and hexagonal arrays and are the same as from the rule of mixture applied in the Vanin model. Other elastic moduli, the transverse Young's modulus * , the axial shear modulus * and the transverse Poisson's ratio * represent the main challenge for the researchers. Here, the predicted values of * , from FE studies are higher than the micromechanical predictions. The current FE model leads also to the relatively high value of the axial shear modulus * that is not predicted by micromechanical models. Lower discrepancy is observed for the transverse shear modulus * . The transverse Poisson's ratio * calculated numerically for the square and hexagonal arrays show higher values than computed according to micromechanical model.

Optimization Problems with the Use of 2D Curves
In optimization problems, distributions of arbitrary design variables can be expressed through n-dimensional curves. This applies not only to the problems of optimizing the shape of the edge of the structure or its central surface but also dimensional optimization (e.g., thickness) or topology/material optimization (distributions: density, the volume fraction of strengthening in the material, the material with different Young's modulus)-see Figure 5. In the latter case, the curves form the boundary between materials with different physical and chemical properties. In cases of simultaneous analysis of optimization, e.g., shape and dimensional, we introduce two types of curves: one describes the variations of the shape of the middle surface, and the other changes the thickness-see Figure 5e. Figure 5. Description of 2D optimization problems with the use of 2D curves: (a) shape optimization-the shape of the middle surface is described by the curve Γ ; (b) shape optimization-the curve describes boundaries of structures Γ = Γ ∪ Γ ∪ Γ ∪ Γ ; (c) dimensional optimization-structure thickness distribution is described by the curve Γ; (d) material optimization-a rectangle is made of three different composite materials (the curves Γ and Γ characterize the division of the area); (e) optimization of the shape of the mid-surface Γ and dimensional optimization (the curve Γ describes the distribution of the thickness).
In this area, optimization of the shape of the external edges of structures (structural shape optimization) is still the most developed in both the 2D and 3D approaches- Figure 5b. Historically, the first work in this field concerned the optimization of the structure's edge-see Ramakrishnan and Francavilla [62]. Currently, the main focus is on finite element (FE) formulations combined with the problems of automatic and even adaptive generation of FE mesh, analytical, semi-analytical or numerical sensitivity analysis (see, e.g., Barthelemy and Haftka [63]) and the development of the effective analytical or numerical methodology of the curve Γ generation.
By the curve, we understand a geometrical figure for which a parametric description can be entered. Generally, it is a set of p points that satisfy the condition p = p(t), where t∈[a, b]. The parametric description of the curve in the n-dimensional space requires the specification of n functions that describe the coordinates of the points of this curve. The method of description Figure 5. Description of 2D optimization problems with the use of 2D curves: (a) shape optimization-the shape of the middle surface is described by the curve Γ; (b) shape optimization-the curve describes boundaries of structures Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 ∪ Γ 4 ; (c) dimensional optimization-structure thickness distribution is described by the curve Γ; (d) material optimization-a rectangle is made of three different composite materials (the curves Γ 1 and Γ 2 characterize the division of the area); (e) optimization of the shape of the mid-surface Γ 2 and dimensional optimization (the curve Γ 1 describes the distribution of the thickness).
In this area, optimization of the shape of the external edges of structures (structural shape optimization) is still the most developed in both the 2D and 3D approaches- Figure 5b. Historically, the first work in this field concerned the optimization of the structure's edge-see Ramakrishnan and Francavilla [62]. Currently, the main focus is on finite element (FE) formulations combined with the problems of automatic and even adaptive generation of FE mesh, analytical, semi-analytical or numerical sensitivity analysis (see, e.g., Barthelemy and Haftka [63]) and the development of the effective analytical or numerical methodology of the curve Γ generation.
By the curve, we understand a geometrical figure for which a parametric description can be entered. Generally, it is a set of p points that satisfy the condition p = p(t), where t ∈ [a, b]. The parametric description of the curve in the n-dimensional space requires the specification of n functions that describe the coordinates of the points of this curve. The method of description (representation) of the curve is directly related to the number of design variables included in the optimization process. Typically, three different variants are used: • Set of nodes' approximates curve Γ (see, e.g., Lee et al. [64]) created by discretization of finite elements, then the number of design variables is large and equal to the number of nodes on the edge Γ, and additionally increases in cases of description concentration of stresses; this method is used very rarely. • Curve Γ is described by elementary analytic functions, in which case it is unambiguously defined by giving a finite fixed number (see Pedersen [65,66]), much smaller than in the previous case; however, there is always a problem with the introduction of the function and assessment, whether it is sufficient to solve a specific optimization problem; this method has unfortunately too little generality. Currently, most of the optimization problems are related to the need for FE discretization of both the curve Γ and the entire structure to determine the value of the objective function. The accuracy of the solution of the problem depends on the accuracy of the approximation of the curve and the numerical model created for the specific task of optimization.
In the further part of the work, we will discuss issues related to two-dimensional curves, i.e., n = 2. We do not intend to discuss here the fundamental aspects of optimization algorithms. They are presented in details in Refs. [67][68][69].

Representation of the Curve Γ through Elementary Functions
The simplest method of building the curve Γ is to express it through analytic functions and a set of parameters. Variations in parameters allow us to optimize the shape. The a priori unknown number of parameters necessary to describe the problem is a distinct disadvantage of such an approach. This methodology can be treated as a starting point for a thorough analysis using spline functions. Already in 1976, Kristiansen and Madsen [70] presented the method of representing the curve in the polar coordinate system in the form: where the function r 0 (θ) defines the basic shape (most often it is an ellipse or a circle-see Figure 6a), and f i (θ) are eigenfunctions that solve the beam vibration equation f IV − λf = 0: where a is a constant resulting from the accepted conditions on the edges θ = 0 and θ = θ 0 , and the parameters t i (real numbers) form a set of design variables.
Materials 2019, 04 , x FOR PEER REVIEW 9 of 24 (representation) of the curve is directly related to the number of design variables included in the optimization process. Typically, three different variants are used: • Set of nodes' approximates curve Γ (see, e.g., Lee et al. [64]) created by discretization of finite elements, then the number of design variables is large and equal to the number of nodes on the edge Γ, and additionally increases in cases of description concentration of stresses; this method is used very rarely. • Curve Γ is described by elementary analytic functions, in which case it is unambiguously defined by giving a finite fixed number (see Pedersen [65,66]), much smaller than in the previous case; however, there is always a problem with the introduction of the function and assessment, whether it is sufficient to solve a specific optimization problem; this method has unfortunately too little generality. Currently, most of the optimization problems are related to the need for FE discretization of both the curve Γ and the entire structure to determine the value of the objective function. The accuracy of the solution of the problem depends on the accuracy of the approximation of the curve and the numerical model created for the specific task of optimization.
In the further part of the work, we will discuss issues related to two-dimensional curves, i.e., n = 2. We do not intend to discuss here the fundamental aspects of optimization algorithms. They are presented in details in Refs. [67][68][69].

Representation of the Curve Γ through Elementary Functions
The simplest method of building the curve Γ is to express it through analytic functions and a set of parameters. Variations in parameters allow us to optimize the shape. The a priori unknown number of parameters necessary to describe the problem is a distinct disadvantage of such an approach. This methodology can be treated as a starting point for a thorough analysis using spline functions. Already in 1976, Kristiansen and Madsen [70] presented the method of representing the curve in the polar coordinate system in the form: where the function ( ) defines the basic shape (most often it is an ellipse or a circle-see Figure  6a), and ( ) are eigenfunctions that solve the beam vibration equation f IV − λf = 0: where a is a constant resulting from the accepted conditions on the edges = 0 and = , and the parameters (real numbers) form a set of design variables. Pedersen et al. [65] extended the number of design variables in the problem of shape optimizing by introducing in Equation (6) unknown parameters of the curve r 0 (θ) defining the basic shape. They proposed a description of this curve as a superellipse whose equation in the Cartesian system has the form: x a n + y b where a, b, n are additional design variables-see Figure 6b. The use of superellipse enables the description of various shapes of the reference curve Γ starting from the straight line (n = 1) and ending with the square (n → ∞). A wider variety of basic shapes is introduced in this way. The curves described by Equations (6) and (8) are used in conjunction with FE models and numerical sensitivity analysis for solving specific optimization tasks.

Generation of Key Point Population
The number of key points and the way they are arranged in a two-dimensional coordinate system (equidistant or not) are directly dependent on the type of optimization problem. It is obvious that the increase in the number of key points increases the size of the task and is unfavourable. On the other hand, increasing the number of key points allows a better representation of the curve Γ. The authors' experiences show that it would be expedient to introduce 5 to 10 basis points to describe one curve. The issue of choosing the location of key points depends primarily on the type of imposed equality or inequality constraints. However, the initialization of the task always starts from the equidistant position of the points. For all key points, we will additionally introduce inequality constraints in the form of: where I specifies the number of key points, and r min and r max are determined depending on the demand; these may be, for example, technological requirements. The values of real numbers r i are determined using a zero-one pseudo-random number generator. We treat the numbers r i as the length of the radii at the end of which the key point P i is located. In the adopted polar coordinate system, it is necessary to define the angular positions of the radii r i . Below are some options for defining them depending on the imposed additional restrictions. It should be emphasized that, in this case, limitations are not boundary conditions imposed at the ends of curve Γ. They are taken into account directly in the methodology of its construction.

No Constraints-Equidistant Key Points
We assume that it is possible to change the angular position of the key points in the range [0 • ,90 • ]. We divide this interval into equally-spaced sections. It is assumed that the first of the generated key points always lie on the OX axis. The vectors start at the origin of the coordinate system. An example of such a generation is shown in Figure 7a. Pedersen et al. [65] extended the number of design variables in the problem of shape optimizing by introducing in Equation (6) unknown parameters of the curve ( ) defining the basic shape. They proposed a description of this curve as a superellipse whose equation in the Cartesian system has the form: where a, b, n are additional design variables-see Figure 6b. The use of superellipse enables the description of various shapes of the reference curve Γ starting from the straight line (n = 1) and ending with the square (n → ∞). A wider variety of basic shapes is introduced in this way. The curves described by Equations (6) and (8) are used in conjunction with FE models and numerical sensitivity analysis for solving specific optimization tasks.

Generation of Key Point Population
The number of key points and the way they are arranged in a two-dimensional coordinate system (equidistant or not) are directly dependent on the type of optimization problem. It is obvious that the increase in the number of key points increases the size of the task and is unfavourable. On the other hand, increasing the number of key points allows a better representation of the curve Γ. The authors' experiences show that it would be expedient to introduce 5 to 10 basis points to describe one curve. The issue of choosing the location of key points depends primarily on the type of imposed equality or inequality constraints. However, the initialization of the task always starts from the equidistant position of the points. For all key points, we will additionally introduce inequality constraints in the form of: where I specifies the number of key points, and and are determined depending on the demand; these may be, for example, technological requirements. The values of real numbers are determined using a zero-one pseudo-random number generator. We treat the numbers as the length of the radii at the end of which the key point Pi is located. In the adopted polar coordinate system, it is necessary to define the angular positions of the radii . Below are some options for defining them depending on the imposed additional restrictions. It should be emphasized that, in this case, limitations are not boundary conditions imposed at the ends of curve Γ. They are taken into account directly in the methodology of its construction.

No Constraints-Equidistant Key Points
We assume that it is possible to change the angular position of the key points in the range [0°,90°]. We divide this interval into equally-spaced sections. It is assumed that the first of the generated key points always lie on the OX axis. The vectors start at the origin of the coordinate system. An example of such a generation is shown in Figure 7a. In the tasks of the curve generation, we are primarily interested in the shape of the geometric figure Γ geometry, and not its location in a particular FE model. In specific numerical cases, it is always possible to translate and rotate the curve to place one (or two) base vectors in the required position.

Generation of Convex Curves
In many optimization problems, the convexity of the generated curve Γ is required. The idea of creating convex curves directly uses the convex polygons' construction algorithm. We randomly generate I + 1 equidistant (in the sense of angular distance) base vectors r i . Next, we create a new sequence of base vectors, starting with the vector with the smallest (largest) length, depending on whether the vector with the smallest (largest) length was generated first in the original sequence. Such a method of selecting the first base vector provides the ability to construct curves for which the base radius on the OX axis may be smaller or larger than the radii on the OY axis. This operation was carried out for the vectors presented in Figure 7a. The new sequence of base vectors (P 0 ,...,P 4 ) is shown in Figure 7b. Next, we check if the vectors form a convex polygon, investigating whether the vector associated with the P i point is above or below the line joining the points (P i−1 , P i+1 ). If it is below the straight line, then we generate a new length so that the end of the vector is beyond the straight line. We start checking the condition of the convexity from the inside, i.e., by examining the position of the vector P 2 relative to the two extreme ones, i.e., P 0 and P 4 . We examine the positions of further points analogously, i.e., dividing subsequent intervals into halves. Of course, the most comfortable is to take an odd number of key points. After completing these operations, repeat the verification of the protuberances relative to the neighbors. In the case shown in Figure 7b, the point P 2 is above the line joining points P 0 and P 4 , but not with respect to the neighboring P 1 and P 3 . Thus, it is necessary to generate a new key point (by generating a new length of the vector) marked in Figure 7b as P 2 Currently, the set of key points forms a convex polygon, and we can construct a convex curve Γ. In computer graphics, the non-uniform rational basis spline (NURBS) mathematical model is commonly used for generating and representing curves and surfaces. In NURBS, control points are entered and the edition of the curve is realized based on the element's control points for Bézier curves or based on spline modeling.

Generation of a convex Curve with Constraints in the Form of a Constant Area Bounded by a Curve Γ
Now, let us discuss the construction of the convex curve Γ with the constraint in the form of the constancy of the A field bounded by this curve and with pre-set boundary conditions in the form of the prescribed values of the angle of tangents at the edges of the curve, i.e., at the key points P 0 and P I . The boundary conditions will be fulfilled as if automatically by determining suitably the positions of the control points necessary for the construction of the curve segments by cubic spline functions. It is a task independent of the generation of key points since such a generation is treated only and exclusively as a problem in the field of geometry. We start the construction of the curve from the generation of the pseudo-random position of the P 0 point on the OX axis, taking into account the Equation (9). Depending on the values of prescribed angles, we distinguish two cases: In the first case (a), the length of the vector r I (the position of the P I point on the OY axis) is determined initially assuming that the searched curve is a quarter of an ellipse with a given area A, i.e.,: The length of this vector should also satisfy Equation (9). The angular positions of the remaining key points are uniquely defined at this stage of the construction, and later they will not be changed anymore. Numerical analysis has shown that the equal angular distances of the base vector positions are not suitable for convex curve construction in the case of a small number of key points. For this reason, we suggest calculating angular positions of base vectors using the information about the arc length of an ellipse with semi-axes (r 0 , r I ). It is assumed that the ellipse segments, obtained after the division of the angle 90 • into the I parts, have equal lengths of arcs. It can be approximately assumed (with sufficient numerical accuracy) that the length of each segment of the ellipse is equal to: In the second case (b), the length of the vector r I (the position of the P I point on the OY axis) is determined by assuming initially that the searched curve is described as the so-called Ferguson curve, which in parametric form is expressed by the formula: where the prime denotes the derivative concerning the parameter t. As previously, it is assumed that the area under the Ferguson curve is equal to A. The angular positions of the key points are determined from a similar condition as for case (a). Next, for the geometrical construction, we generate pseudo-randomly the remaining lengths of base vectors, demanding that two conditions be satisfied: 1.
Convexity of a polygon stretched on basis vectors, 2.
Boundary conditions at the ends of the curve.
The numerical implementation of the first condition is analogous to that described in the previous chapter. The second condition imposes only geometric constraints on the area of possible positions of base vectors, as shown in Figure 8. The length of this vector should also satisfy Equation (9). The angular positions of the remaining key points are uniquely defined at this stage of the construction, and later they will not be changed anymore. Numerical analysis has shown that the equal angular distances of the base vector positions are not suitable for convex curve construction in the case of a small number of key points. For this reason, we suggest calculating angular positions of base vectors using the information about the arc length of an ellipse with semi-axes ( , ). It is assumed that the ellipse segments, obtained after the division of the angle 90° into the I parts, have equal lengths of arcs. It can be approximately assumed (with sufficient numerical accuracy) that the length of each segment of the ellipse is equal to: In the second case (b), the length of the vector (the position of the PI point on the OY axis) is determined by assuming initially that the searched curve is described as the so-called Ferguson curve, which in parametric form is expressed by the formula: where the prime denotes the derivative concerning the parameter t. As previously, it is assumed that the area under the Ferguson curve is equal to A. The angular positions of the key points are determined from a similar condition as for case (a). Next, for the geometrical construction, we generate pseudo-randomly the remaining lengths of base vectors, demanding that two conditions be satisfied: 1. Convexity of a polygon stretched on basis vectors, 2. Boundary conditions at the ends of the curve.
The numerical implementation of the first condition is analogous to that described in the previous chapter. The second condition imposes only geometric constraints on the area of possible positions of base vectors, as shown in Figure 8.  After selecting the positions of base vectors, it is necessary to check the condition of the constant area bounded by a polygon. The relationship for calculating the triangle area OP i−1 P i is used here: If the sum of the fields of individual triangles is not equal to A, then we reduce or increase the value of the radius r I and repeat the operations of selecting the positions of the vectors starting from the correction of the angular positions of the base vectors. The construction of the curve that meets the condition of the convexity and constancy of the area is a set of geometrical observations based on the idea of rejecting solutions (positions) that do not meet the imposed conditions (unsuccessful trials). This method is simple for numerical algorithms and, in our opinion, has large generalisations in the sense of being easily adaptable to a range of similar problems.

2D Homogenization Problem
The homogenization of 2D microcomposites has been presented for several examples of RVE, i.e., a structure with short fibre (Figure 9a), a bilayered structure (Figure 9b), a structure with a circular reinforcement (Figure 9c), and a structure with a triangular shape of reinforcement (Figure 9d). The results of numerical homogenization and values calculated from micromechanical models have been presented and discussed. The problem of the inclusion shape and its influence on the effective properties of heterogeneous media is especially emphasized in the literature, e.g., in the paper by Savvas et al. [71]. the idea of rejecting solutions (positions) that do not meet the imposed conditions (unsuccessful trials). This method is simple for numerical algorithms and, in our opinion, has large generalisations in the sense of being easily adaptable to a range of similar problems.

2D Homogenization Problem
The homogenization of 2D microcomposites has been presented for several examples of RVE, i.e., a structure with short fibre (Figure 9a), a bilayered structure (Figure 9b), a structure with a circular reinforcement (Figure 9c), and a structure with a triangular shape of reinforcement ( Figure  9d). The results of numerical homogenization and values calculated from micromechanical models have been presented and discussed. The problem of the inclusion shape and its influence on the effective properties of heterogeneous media is especially emphasized in the literature, e.g., in the paper by Savvas et al. [71]. The numerical homogenization procedure was applied to analyze 2D composites with the uniform distribution of short fibres (Figure 9a), with the bilayered structure (Figure 9b), a structure with a circular reinforcement (Figure 9c), and with the triangular shape of reinforcement (Figure 9d). The assumed material data are listed in Table 5, and volume fraction of reinforcement for the short fibres, circular reinforcement and the bilayered structure was v = 32%, whereas, for the triangular shape of the reinforcement, v = 10%. The analytical solution for the problem of 2D composites with the uniform distribution of short fibres is presented, e.g., in the paper by Fuji and Zako [72]. Here, a short fibre having ratio d/l = 0.5 was analyzed (Figure 9a). The uni-axial loading in the elastic regime was applied for all 2D RVEs- Figure 9 and plane strain elements (CPE4 in ABAQUS package) were used.  The numerical homogenization procedure was applied to analyze 2D composites with the uniform distribution of short fibres (Figure 9a), with the bilayered structure (Figure 9b), a structure with a circular reinforcement (Figure 9c), and with the triangular shape of reinforcement (Figure 9d). The assumed material data are listed in Table 5, and volume fraction of reinforcement for the short fibres, circular reinforcement and the bilayered structure was v f = 32%, whereas, for the triangular shape of the reinforcement, v f = 10%. The analytical solution for the problem of 2D composites with the uniform distribution of short fibres is presented, e.g., in the paper by Fuji and Zako [72]. Here, a short fibre having ratio d/l = 0.5 was analyzed (Figure 9a). The uni-axial loading in the elastic regime was applied for all 2D RVEs- Figure 9 and plane strain elements (CPE4 in ABAQUS package) were used. In the 2D analysis of microcomposites, the effective Young's modulus E * in the loading direction as the ratio of average stresses σ ij to average strains ε ij (Equations (2) and (3)) was computed numerically. The E * evaluated numerically was compared with rule of the mixture (ROM), Halpin-Tsai (H-T) model for short fibers, Cox model [73] and micromechanical model presented by Strzelecki [74] (see Equations (A2)-(A5)). The results are listed in Table 6. For all analyzed cases, the effective Young's modulus according to the rule of mixture is higher than the numerical values. The numerical homogenization results are close to the Halpin-Tsai model for short fibres, whereas the Cox model gives lower values of E * . The FEM values of E * for the bilayered stucture are very close to the micromechanical models. Similarly, as in cases (a) and (c)- Figure 9, the numerical results of RVE with the triangular shape of the reinforcement are much lower than the values from the rule of mixture. The applied micromechanical models give higher results, but the difference is not so high. However, it should be noted that the shape of the reinforcement is very complicated compared to those adopted in micromechanical models. This evokes other distributions of stresses and strains in the representative cell and thus causes differences in the relation to simple micromechanical models. Practically, such a complicated shape excludes the possibility of using these models due to difficulties in determining the diameter of the reinforcement. Here, the reporting results in Table 6 assumed d/l = 0.075 in analytical computations. Table 6. Effective Young's modulus of microcomposites (material data in Table 5). The values of the equivalent Young's modulus obtained from the numerical analysis for the 2D problems presented above were significantly different from the results calculated based on the rule of mixture. The application of the simple rule of mixture overestimates the value of the effective Young's modulus (see Table 6). The use of micromechanical models is appropriate in the case of simple shape reinforcement, i.e., an easy determination of the reinforcement diameter and length. In the case of a complicated shape of the reinforcement (see Figure 9d), phenomenological models should not be used because the determination of d/l ratio is problematic.

Isoperimetric Problem-Verification of the Accuracy of Numerical Solutions
If we impose equality conditions limiting the sought functions, then the optimization of the shape (looking for the extremum of a function dependent on the function defining the shape of the structure's edge) with imposed constrained on the shape function (in the implicit and functional form) are the classic, isoperational variational problems. The solutions of these tasks are sought in the class of differentiable functions with a given metric and satisfying the imposed equality constraints-see Tatarkiewicz [75]. One of the most well-known issues in this field is the following problem: Find a closed plane curve of a given perimeter L (refers to the length of the boundary connecting two points A and B) which encloses the maximal Area.
The solution is a circle (the class of curves is limited to smooth curves) with an unknown location of the center of the circle and radius. These quantities are determined from the boundary conditions.
The above problem is very convenient to check to what extent the application of Bezier's spline functions is sufficient to approximate continuous convex functions. However, to obtain a specific numerical solution, it is necessary to specify constants. We assume that the center of the circle is located at the beginning of the Cartesian coordinate system and the arc length of the circle L connecting the points A and B is πr, where r is half the length of the segment AB. The solution to the problem is, therefore, a semi-circle with an area of πr 2 /2. Due to the symmetry of the problem concerning the y-axis, we only look for half of the curve. Finally, we formulate the optimization problem in the following way: where Area means the field under the convex curve Γ with the arc length L, and η the penalty parameter. This amount was assumed in the calculations to be equal to 250. It was also assumed that the curves sought are always convex. Examples of a generation of nine different populations are shown in Figure 10. The area of the surface under each of the generalized curves is equal to the area of the circle. The solution is a circle (the class of curves is limited to smooth curves) with an unknown location of the center of the circle and radius. These quantities are determined from the boundary conditions.
The above problem is very convenient to check to what extent the application of Bezier's spline functions is sufficient to approximate continuous convex functions. However, to obtain a specific numerical solution, it is necessary to specify constants. We assume that the center of the circle is located at the beginning of the Cartesian coordinate system and the arc length of the circle L connecting the points A and B is πr, where r is half the length of the segment AB. The solution to the problem is, therefore, a semi-circle with an area of πr 2 /2. Due to the symmetry of the problem concerning the y-axis, we only look for half of the curve. Finally, we formulate the optimization problem in the following way: Find where Area means the field under the convex curve Γ with the arc length , and η the penalty parameter. This amount was assumed in the calculations to be equal to 250. It was also assumed that the curves sought are always convex. Examples of a generation of nine different populations are shown in Figure 10. The area of the surface under each of the generalized curves is equal to the area of the circle. The results of calculations for a population of 50 individuals are shown in Figure 11. In the calculations, it has been assumed that functions must always pass through a point with coordinates (4.5,0.0). Furthermore, the possible coordinate variations in the vertical direction were limited, i.e., 0.0 ≤ y ≤ 1.2 × 4.5 = 5.4.
As can be seen, in Figure 11 and Table 7, there is a good agreement of theoretical and numerical solutions, both in the case of using genetic algorithms (GA) as well as a modified evolutionary strategy (MEA). The results of calculations for a population of 50 individuals are shown in Figure 11. In the calculations, it has been assumed that functions must always pass through a point with coordinates   Figure 11. Optimal curves-GA and MEA (black color) and obtained using the optimization procedure (red color).
The number of generations necessary to obtain optimal solutions is small-see Figure 11. It is mainly due to the imposition of a constraint in the form of a demand to fulfil the condition that curves Γ are to be always convex. In this way, generated curves that do not meet this condition are not taken into account in the results presented in Figure 12. As can be seen, in Figure 11 and Table 7, there is a good agreement of theoretical and numerical solutions, both in the case of using genetic algorithms (GA) as well as a modified evolutionary strategy (MEA). The number of generations necessary to obtain optimal solutions is small-see Figure 11. It is mainly due to the imposition of a constraint in the form of a demand to fulfil the condition that curves Γ are to be always convex. In this way, generated curves that do not meet this condition are not taken into account in the results presented in Figure 12. The number of generations necessary to obtain optimal solutions is small-see Figure 11. It is mainly due to the imposition of a constraint in the form of a demand to fulfil the condition that curves Γ are to be always convex. In this way, generated curves that do not meet this condition are not taken into account in the results presented in Figure 12.

Shape of Fibre Bundles
Unidirectional composites are not reinforced with a single fiber but a bundle of fibers embedded in the matrix. The shape of the cross section of the fiber bundle determined as the edge of the area occupied by the fiber Γ = Ω can be arbitrary, i.e., it can form a circle, an ellipse, etc. Using

Shape of Fibre Bundles
Unidirectional composites are not reinforced with a single fiber but a bundle of fibers embedded in the matrix. The shape of the cross section of the fiber bundle determined as the edge of the area occupied by the fiber Γ = ∂Ω f can be arbitrary, i.e., it can form a circle, an ellipse, etc. Using the numerical homogenization (see Section 2 and [31]), we can determine the mechanical properties of an elementary cell made of fiber and matrix. The area occupied by the cell is denoted as Ω = Ω f ∪ Ω m . The shape of the edge of the area occupied by the fiber has a significant effect on the homogenized mechanical properties for the elementary cell. Thus, it also affects the stiffness of the laminate. The purpose of the optimization is to determine the edge shape Γ = ∂Ω f , to minimize the function of the total internal energy of the RVE in the Lagrange form, i.e.,: with a constraint that the area enclosed by the curve Γ is constant, i.e.,: where g Γ defines a function describing the curve. Equation (15) can be written in the following explicit form: We approximate the curve Γ using the spline functions by entering five key points. The condition in the form of Equation (16) will be written in the dimensionless form during the generation of key points. In the numerical analysis, we consider only 1/4 of the elementary cell. Therefore, additional restrictions on the derivative value at the ends of the curve were imposed.
The elemental cell was modeled using triangular FE type NKTP 1 (plane stress state-NISA package). On the edge of the elementary cell, external exertion was applied in the form of a unit displacement acting along the x-axis. A circle was assumed as the original shape of the fiber bundle. The results of solving Equations (15) and (16) are shown in Figure 13. In the numerical computations, the algorithm of the modified evolutionary strategy (MEA) was used. The shape of the fiber bundle depends on the conditions limiting the form of the curve Γ understood in the sense of its convexity or not.
where defines a function describing the curve. Equation (15) can be written in the following explicit form: We approximate the curve Γ using the spline functions by entering five key points. The condition in the form of Equation (16) will be written in the dimensionless form during the generation of key points. In the numerical analysis, we consider only 1/4 of the elementary cell. Therefore, additional restrictions on the derivative value at the ends of the curve were imposed.
The elemental cell was modeled using triangular FE type NKTP 1 (plane stress state-NISA package). On the edge of the elementary cell, external exertion was applied in the form of a unit displacement acting along the x-axis. A circle was assumed as the original shape of the fiber bundle. The results of solving Equations (15) and (16) are shown in Figure 13. In the numerical computations, the algorithm of the modified evolutionary strategy (MEA) was used. The shape of the fiber bundle depends on the conditions limiting the form of the curve Γ understood in the sense of its convexity or not. Figure 13. Planar representative volume element-initial (circular) and optimal shapes of fibre bundles.
The presented numerical solution draws new possibilities to optimize the mechanical properties of the material by forming a bundle of reinforcement fibers. It is possible from a technological point of view. However, attention should be paid to the issue of the interfacial surface, which was not included in the considerations. The shape of the fiber bundle may have some influence on the initialisation process of the micro failure at the fiber/matrix boundary and may be responsible for the formation of macrocracks and destruction of the material. The presented numerical solution draws new possibilities to optimize the mechanical properties of the material by forming a bundle of reinforcement fibers. It is possible from a technological point of view. However, attention should be paid to the issue of the interfacial surface, which was not included in the considerations. The shape of the fiber bundle may have some influence on the initialisation process of the micro failure at the fiber/matrix boundary and may be responsible for the formation of macrocracks and destruction of the material.

Distribution of the Reinforcement
The functionally gradient materials (FGM) distribute the material functions throughout the material body to achieve the maximum mechanical properties. FGMs are characterized as structures having different properties in the thickness z-direction in the form proposed by Qian et al. [39]: where p denotes a generic material property like modulus, p t and p b denote the corresponding properties of the top and bottom faces of the plate, respectively, and V(z) is a function describing the distribution of reinforcement density in the thickness direction. For FGM, the function V(z) is defined as: Other definitions of the V(z) can be found, e.g., in reference [40]. Other propositions consider the use of Equations (18) and (19) for degradation in a different direction than thickness. First density grading in the width direction is as follows: where ρ l and ρ t are the mass densities of the leading edge and trailing edge, respectively; w denotes the width of the panel. Second density grading in the length direction is as follows: where ρ r and ρ t are the mass densities of the root edge and tip edge, respectively; L denotes the length of the panel. Next definition of density diagonal grading across the xy plane (in two directions) is as follows: where ρ u and ρ l are the mass densities of upper and lower bounds, respectively. A model for characterizing the mechanical properties of functionally graded materials (FGMs) with the regular polygonal cross-section is also developed by introducing the power-law rule-see Ref. [41]. However, modelling of their mechanical behaviour still leads to many problems, particularly in the description of 1D or 2D structures, such as beams, plates or shells. The broader discussion of those problems and solutions can be found in Muc et al. [76][77][78][79][80]. Let us consider a composite beam reinforced by short fibres, loaded at the center ( Figure 14three-point bending test) and having a variable density of fibres ρ f (x) along the length. The present methodology is similar to the FGM description. In the current problem, the density grading in the length direction is presented. The objective of the optimisation is as follows: to find the minimal mass of the beam satisfying the equality constraint imposed on the displacement parameter and inequality (upper and lower bounds) constraint-on the fibre volume density fraction Vf-the strict formulation is given by Banichuk et al. [81]. The optimisation problem deals with the topology optimisation since we are looking for the optimal material distribution. The results are presented in Figure 15 and show a good correlation between numerical and analytical studies. Upper values (red lines) in Figure 15 were calculated for the lower value of the assumed beam deflection. The objective of the optimisation is as follows: to find the minimal mass of the beam satisfying the equality constraint imposed on the displacement parameter and inequality (upper and lower bounds) constraint-on the fibre volume density fraction V f -the strict formulation is given by Banichuk et al. [81]. The optimisation problem deals with the topology optimisation since we are looking for the optimal material distribution. The results are presented in Figure 15 and show a good correlation between numerical and analytical studies. Upper values (red lines) in Figure 15 were calculated for the lower value of the assumed beam deflection.
The objective of the optimisation is as follows: to find the minimal mass of the beam satisfying the equality constraint imposed on the displacement parameter and inequality (upper and lower bounds) constraint-on the fibre volume density fraction Vf-the strict formulation is given by Banichuk et al. [81]. The optimisation problem deals with the topology optimisation since we are looking for the optimal material distribution. The results are presented in Figure 15 and show a good correlation between numerical and analytical studies. Upper values (red lines) in Figure 15 were calculated for the lower value of the assumed beam deflection.

Conclusions
The present paper discussed the application of numerical homogenization and optimization in the studies of micro-and nanocomposites. The issue of boundary conditions, reinforcement shape and distribution, and form of representative volume element is presented taking into account also the topic of 2D curves optimization.
The conducted analyses and comparisons are summarized as follows: • The method of description (representation) of the 2D curve related to the number of design variables included in the optimization process was presented.

•
There was a good agreement of theoretical and numerical solutions, both in the case of using genetic algorithms (GA) as well as a modified evolutionary strategy (MEA) for the shape optimization of the single fibre.

Conclusions
The present paper discussed the application of numerical homogenization and optimization in the studies of micro-and nanocomposites. The issue of boundary conditions, reinforcement shape and distribution, and form of representative volume element is presented taking into account also the topic of 2D curves optimization.
The conducted analyses and comparisons are summarized as follows: • The method of description (representation) of the 2D curve related to the number of design variables included in the optimization process was presented.

•
There was a good agreement of theoretical and numerical solutions, both in the case of using genetic algorithms (GA) as well as a modified evolutionary strategy (MEA) for the shape optimization of the single fibre.

•
In shape optimization of fibre bundles, a circle was assumed as the original shape. In the numerical computations, the algorithm of modified evolutionary strategy (MEA) was used. The shape of the fibre bundle depends on the conditions limiting the form of the curve Γ understood in the sense of its convexity or not.

•
The composite beam reinforced by short fibres, loaded at the centre and having a variable density of fibres along the length, was considered, being an example of FGM. The optimization problem deals with topology optimization since we were looking for the optimal material distribution. The results showed a good correlation between numerical and analytical studies.

•
Some examples of numerical homogenization of various 2D and 3D RVE have been studied to present the influence of boundary conditions, form of RVE, shape and distribution of the reinforcement on the effective material properties.

•
The comparison of numerical homogenization and micromechanical models showed that micromechanical models are appropriate in the case of the simple shapes of reinforcement. • However, for complicated shapes of the reinforcement, the application of the numerical homogenization is crucial because the determination of characteristic relations such as d/l ratio is problematic. Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest.