Introduction to Macroscopic Optimal Design in the Mechanics of Composite Materials and Structures

: The main goal of building composite materials and structures is to provide appropriate a priori controlled physico-chemical properties. For this purpose, a strengthening is introduced that can bear loads higher than those borne by isotropic materials, improve creep resistance, etc. Composite materials can be designed in a different fashion to meet speciﬁc properties requirements.Nevertheless, it is necessary to be careful about the orientation, placement and sizes of different types of reinforcement. These issues should be solved by optimization, which, however, requires the construction of appropriate models. In the present paper we intend to discuss formulations of kinematic and constitutive relations and the possible application of homogenization methods. Then, 2D relations for multilayered composite plates and cylindrical shells are derived with the use of the Euler–Lagrange equations, through the application of the symbolic package Mathematica. The introduced form of the First-Ply-Failure criteria demonstrates the non-uniqueness in solutions and complications in searching for the global macroscopic optimal solutions. The information presented to readers is enriched by adding selected review papers, surveys and monographs in the area of composite structures.


Introduction-General Characteristics and Applications of Composite Materials (CM)
The idea of composite materials arose from the need to combine different materials in order to overcome the shortcomings of each one of them. In this way, the properties of the composites are better than the properties of the phases that form them. One of the phases (materials) is called the matrix, and the other, reinforcing the composite, is called reinforcement.
The progress that has been made in the development of CM and its impact on civilization transformations is fast and unusual [1]. It mainly results from the numerous advantages of these materials, the most important of which are: (a) low specific weight ρ, (b) high resistance of CM with a polymer or ceramic matrix for corrosion, (c) low thermal conductivity; this value is much lower than for steel, (d) ceramic CM resistance to temperature; their mechanical properties do not change even at 2000 • C, (e) good CM resistance with a polymer matrix for cracking; using the critical stress factor K c as a measure, this resistance is comparable to the resistance of medium carbon steels and Al alloys, and when evaluated by comparing the critical rate of G c release rate, the fracture toughness corresponds to that of soft steels or Ti alloys; in the case of ceramic matrix composites, the above coefficients are much lower, (f) high CM resistance with polymer matrix for local damage; steel structures are very sensitive to local damage, e.g., poor weld performance, (g) technological ease of manufacturing objects (constructions) with complex shapes in single operations, (h) composites are in principle a composition of all types of materials, e.g., a combination of materials, one of which has piezo-electric properties; the other, electro-luminescent, provides a composite with piezo-luminescent properties, shining under the action of an applied force.

•
The physical relations can be built separatelyas a description of the state of deformation (Section 3), taking also into account homogenization problems (Section 4).

•
The formulation of the relation strain-deformation (displacements) is completely independent on the physical relations.

•
The optimization problems can be formulated and solved at the macroscale level of formulations, and the sense of approximations used in the physical relations and the strain displacement relations can be verified only by the comparison with experimental results or other numerical solutions.

Material Selection vs. Material Design-The Necessity of Optimization
Up until the end of the 20th century, human needs have been met by designing and producing constructions for which different materials have been selected. In this way, the basic emphasis was placed on the optimal shaping of constructional features, such as shape, dimensions and weight, assuming a priori that material properties were strictly defined or could be selected within a strictly defined scope. This tendency is clearly outlined in Ashby's monograph [2].
Nowadays, however, the development of materials science and materials engineering indicates a different direction of material engineering activities (see, e.g., Rühle et al. [3]), called material design. This term means the adaptation of materials from the chemical composition, component phases and nano-microstructure to the properties required for their application. In the not too distant future, traditional, empirical methods of introducing new materials will be more and more supplemented with theoretical predictions.
Differences in the concepts of material design and material selection are shown in Figure 1. In the case of material selection, the construction plays a superior role in relation to the material (arrow direction), because the construction is based on an existing material or group of materials. In the second case (material design), the situation is reversed. The overriding issue is the material that we design properly from scratch. The relationship with the structure is that we demand the fulfillment of specific functional properties by the material. However, they are specific only for a single structure. For another construction, we can design a completely different material.
to the material (arrow direction), because the construction is based on an existing material or group of materials. In the second case (material design), the situation is reversed. The overriding issue is the material that we design properly from scratch. The relationship with the structure is that we demand the fulfillment of specific functional properties by the material. However, they are specific only for a single structure. For another construction, we can design a completely different material.

Figure 1. Concepts of materials selection and design.
Usually, we intend to allow the pressure vessel to carry high pressures, and we also require it to be light (easy to transport). In the era of the existence of isotropic materials, the constructors encountered a barrier of material properties. For a given geometry of the pressure vessel (length, diameter), the value of the pressure transferred and the thickness of the vessel were related to the boiler codes with allowable strength (material property). For the highest possible strength (limited group of materials), the increase in pressure is always associated with the increase in the thickness of the tank (increase in weight). In the material design situation, it is possible to demand the increase of the allowable strength by designing a completely new fibrous CM (designing yet non-existent fibers, special finishing and matrix) and in this way increase the transferred pressure with a low weight of the structure. Currently, thanks to CM, very light pressure vessels are created for pressures of several hundred atmospheres. However, it is necessary to clearly distinguish the material design for specific structures from the continuously growing set of materials available for selection. Ashby [2] gives the example of a vacuum cleaner whose weight and price have been significantly reduced compared to the initial situation with the aid of the discovery of polymers.
In the case of material design, its structure and properties (physicochemical, nuclear, electrical, mechanical, technological, etc.) have a great influence not only on the type of structure (beam, plate, coating, load method, boundary conditions), but also on the technological process ( Figure 1).

Material Selection
Material Design Usually, we intend to allow the pressure vessel to carry high pressures, and we also require it to be light (easy to transport). In the era of the existence of isotropic materials, the constructors encountered a barrier of material properties. For a given geometry of the pressure vessel (length, diameter), the value of the pressure transferred and the thickness of the vessel were related to the boiler codes with allowable strength (material property). For the highest possible strength (limited group of materials), the increase in pressure is always associated with the increase in the thickness of the tank (increase in weight). In the material design situation, it is possible to demand the increase of the allowable strength by designing a completely new fibrous CM (designing yet non-existent fibers, special finishing and matrix) and in this way increase the transferred pressure with a low weight of the structure. Currently, thanks to CM, very light pressure vessels are created for pressures of several hundred atmospheres. However, it is necessary to clearly distinguish the material design for specific structures from the continuously growing set of materials available for selection. Ashby [2] gives the example of a vacuum cleaner whose weight and price have been significantly reduced compared to the initial situation with the aid of the discovery of polymers.
In the case of material design, its structure and properties (physicochemical, nuclear, electrical, mechanical, technological, etc.) have a great influence not only on the type of structure (beam, plate, coating, load method, boundary conditions), but also on the technological process (Figure 1). Figure 1 demonstrates that in order to design materials for a particular structure, it is necessary to analyze the mutual dependencies existing between the four main components: construction-technology-structure-material.
Modern methods of design or material selection are using optimization methods to an increasing extent. This is due, on the one hand, to the increasing number of components (criteria) necessary to take into account to ensure the a priori imposed general and specific requirements in the project under development and, on the other hand, to the development and application possibilities of modern computer-aided design methods. The implementation of these concepts also triggers a greater attention to the issues of optimal design of technological processes used to produce a material that fulfills specific functions in a specific construction (e.g., load transfer, deformation, stiffness, etc.). Usually, multi-criteria optimization methods should be used in the selection or design of materials, first of all in order to maintain an appropriate compromise between the components presented in Figure 1.
In the case of modern materials, including composite materials, the problem of material design is associated with the need to go down to a lower level of description of interactions, understood as the transition from the macro state to the meso and micro-tonano-structures (scale of atomic interactions) in order to obtain appropriate (consistent with the results of experiments) structural responses of the system on a macro scale. In the case of fibrous composite materials, experimental investigations significantly exceed the theoretical modeling of their properties, mainly due to the complexity of processes occurring at the interface between the fiber and the matrix (local non-linear effects). Due to the complexity of the physical and chemical effects, the development of microcracks and the variability of the reinforcement cross-sections, we notice in the experimental evaluation of the mechanical properties a significant scatter of the results in relation to the average value. This leads directly to the need to consider not only purely deterministic approaches in the description of mechanical properties and the optimization of structures made of composite materials, but also the technological processes of their implementation. The above difficulties lead to the use, in relation to the description and analysis of mechanical properties of structures made of fibrous CM (microcomposites), of the homogenization method. There is a rich, long-term documentation in the global literature demonstrating the advisability of applying such an assumption in both scientific and engineering practice.
Currently, as in the past, the process of optimal design is an interdisciplinary activity and can still be described by the triangle proposed by Eschenauer [4] (Figure 2). However, it should be clearly emphasized that the problem of the optimization of composite structures and their production processes is an area even younger than the mechanics of composites. In addition, there are many new problems in this area of optimization, completely different from the tasks of optimizing structures made of isotropic materials.  Due to the above facts, the majority of problems to be solved concern fundamental issues, i.e.: -modeling the optimization problem by taking into account (usually the simplest) models of failure of structures and models of technological processes, -attempts to take into account statistical distributions of experimental results-statistical approach, anti-optimization, fuzzy set theory, -selection of appropriate design variables and definitions of design spaceespecially appropriate for introducing various constraints in the form of equality or inequality, -choosing the right optimization algorithm, which becomes extremely important in the case of many design variables, the existence of many local extremes of the objective function and the lack of convexity for many optimization problems, Due to the above facts, the majority of problems to be solved concern fundamental issues, i.e.,: modeling the optimization problem by taking into account (usually the simplest) models of failure of structures and models of technological processes, -attempts to take into account statistical distributions of experimental results-statistical approach, anti-optimization, fuzzy set theory, -selection of appropriate design variables and definitions of design spaceespecially appropriate for introducing various constraints in the form of equality or inequality, -choosing the right optimization algorithm, which becomes extremely important in the case of many design variables, the existence of many local extremes of the objective function and the lack of convexity for many optimization problems, -the necessity to search for solutions based on numerical methods (e.g., finite element method, finite volume method), which results mainly from the lack of analytical solutions for optimization tasks in the field of composite mechanics-see the quotation from Tsai, Hahn [5]: It is unfortunate that the use of composite materials is limited or penalized by the nonavailability of analytical tools. It is important to understand how anisotropy and nonhomogenity arise in composite laminates and to what degree they can be manipulated to perform functions not possible with conventional materials.
The concept of topological optimization includes a full set of design variables that combine dimension, shape, topological and material-related design variables. However, the above generalization should not be confused with a strictly understood definition of topology (see, e.g., [6,7]), but understood strictly in light of the theory of optimization of topological variables.
In the optimization of composite materials or constructions made of composite materials, the layout optimization of composite structures is the most adequate, instead of topological optimization. This is due to the fact that topological optimization is mainly associated only with the term of material homogenization (see [6]). In my opinion, the most reasonable approach seems to be the use of the definition: the optimization of the structure of composite materials with the simultaneous distinction of the level of analysis (the length scale), i.e.,it is possible to optimize at a specific level (see Figure 3 This is necessary because on each of the above-mentioned levels of analysis CM introduce different types of design variables, including both the geometry and the type of reinforcement, matrix and interfacial layer or the orientation of individual layers, the individual thickness of individual layers and their material properties. The set of design variables describing the structure of CM is an arrangement of geometrical and physical variables that characterize (at a given level of analysis) unambiguously the physical properties of CM. The next difference between the definitions proposed in this work and the Bendsøe approach [6,7] lies in the possibility of describing the structure by other design variables than those that characterize the problem's physics, namely dimensional and shape optimization variables describing the geometry of the problem. In this sense, we use two types of design variables: • physical (material), representing the CM from which the structure is made; • geometric, characterizing the geometry of the structure. This is necessary because on each of the above-mentioned levels of analysis CM introduce different types of design variables, including both the geometry and the type of reinforcement, matrix and interfacial layer or the orientation of individual layers, the individual thickness of individual layers and their material properties. The set of design variables describing the structure of CM is an arrangement of geometrical and physical variables that characterize (at a given level of analysis) unambiguously the physical properties of CM. The next difference between the definitions proposed in this work and the Bendsøe approach [6,7] lies in the possibility of describing the structure by other design variables than those that characterize the problem's physics, namely dimensional and shape optimization variables describing the geometry of the problem. In this sense, we use two types of design variables:  physical (material), representing the CM from which the structure is made;  geometric, characterizing the geometry of the structure.

Basic Definitions and Notations
The problems discussed in the paper concern micro-composites, i.e., composite materials with a reinforcement embedded in a polymeric resin.
The structural analysis of elements composed of composite materials is based on the theory of anisotropic elasticity. For a linear anisotropic material, the generalized Hooke's law relates the stress σi and the strain tensor εj [5]: σi = Cijεj, i,j = 1, 2, ..., 6 (1) or, in the inverse form: εi = Sijσj, i,j = 1, 2, ..., 6 (2) where: Cij are components of the stiffness matrix (Sij-the compliance or flexibility matrix). The relations are written in the simplified, contracted single-subscript notation for the stress and strain components and a double-subscript notation for the elastic con-

Basic Definitions and Notations
The problems discussed in the paper concern micro-composites, i.e., composite materials with a reinforcement embedded in a polymeric resin.
The structural analysis of elements composed of composite materials is based on the theory of anisotropic elasticity. For a linear anisotropic material, the generalized Hooke's law relates the stress σ i and the strain tensor ε j [5]: σ i = C ij ε j , i,j = 1, 2, ..., 6 (1) or, in the inverse form: ε i = S ij σ j , i,j = 1, 2, ..., 6 where: C ij are components of the stiffness matrix (S ij -the compliance or flexibility matrix). The relations are written in the simplified, contracted single-subscript notation for the stress and strain components and a double-subscript notation for the elastic constraints (so-called Kelvin-Voigt's notation) [5]. They are introduced in an orthogonal system of coordinates ( Figure 4). 23 , ε 5 = 2ε 13 , ε 6 = 2ε 12 (3) 7 of 25 system of coordinates ( Figure 4). , , ,  Due to the assumed symmetry of the stress and the strain tensors, the number of the independent non-zero material coefficients in the stiffness matrix Cij is equal to 36. A further reduction in the number of their values is possible if additional physical and material assumptions are introduced into our analysis. The values of the material coefficients depend on the reference coordinate system, since in the general tensor notation the material coefficient matrix is a fourth-order tensor, and, therefore, their coefficients follow tensor transformation rules. Now, let us assume that the constitutive relations (1) are formulated in the local (material) Cartesian coordinate system.
For composite materials having an elastic potential function and three symmetry planes that are mutually orthogonal in the material system of coordinates, the stiffness matrix has the following structure-Vinson Such a material is called an orthotropic material (the shortening of orthogonal anisotropy). In this case the number of independent material parameters is reduced to nine engineering constants (E1, E2, E3, G12, G13, G23, ν12, ν13, ν23). The above form of the stiffness matrix can be easily adopted not only for the analysis of laminated multilayered materi- Due to the assumed symmetry of the stress and the strain tensors, the number of the independent non-zero material coefficients in the stiffness matrix C ij is equal to 36. A further reduction in the number of their values is possible if additional physical and material assumptions are introduced into our analysis. The values of the material coefficients depend on the reference coordinate system, since in the general tensor notation the material coefficient matrix is a fourth-order tensor, and, therefore, their coefficients follow tensor transformation rules. Now, let us assume that the constitutive relations (1) are formulated in the local (material) Cartesian coordinate system.
For composite materials having an elastic potential function and three symmetry planes that are mutually orthogonal in the material system of coordinates, the stiffness matrix has the following structure-Vinson [8]: where: Such a material is called an orthotropic material (the shortening of orthogonal anisotropy). In this case the number of independent material parameters is reduced to nine engineering constants (E 1 , E 2 , E 3 , G 12 , G 13 , G 23 , ν 12 , ν 13 , ν 23 ). The above form of the stiffness matrix can be easily adopted not only for the analysis of laminated multilayered materials but also for the investigation of structures made of porous functionally-graded materials [9] or nanostructures embedded in a matrix [10]. The identical formulation can also be easily extended to the problems of structural deformations of constructions with piezoelectric sensors/actuators [11,12]. However, in the latter case, the stiffness matrix (1) or (4) is also extended by adding the Maxwell relations.

Homogenization of Mechanical Properties
Polymer composites reinforced with unidirectional fibers (1D composites) or textile composites (2D or 3D fiber structures) are heterogeneous (inhomogeneous) materials, since they are made of two or more constituents with different physicochemical properties. A microscopically inhomogeneous composite material can be idealized as a macroscopically homogeneous continuum when the behavior of engineering structures made of the material can be satisfactorily retained. Such idealization can be realized over a representative sample of the composite material. The selection of the dimensions of a representative volume is imperative. The representative volume must be sufficiently large compared to the scale of the microstructure, so that it contains a sufficient number of individual constituents (see, e.g., Neto et al. [13,14]). In addition, it should be entirely typical of the whole composite structure on average (see, e.g., Figure 5). It is assumed that the representative volume contains all microstructural features, i.e.,: • a set of mechanical parameters describing the mechanical properties of constituents, • geometrical parameters characterizing the dimensions of a reinforcement and its location in a representative cell, • a position of a representative cell in a global coordinate system; it may be expressed by Euler's angles in three dimensions or by an angle θ in two dimensions.
since they are made of two or more constituents with different physicochemical properties. A microscopically inhomogeneous composite material can be idealized as a macroscopically homogeneous continuum when the behavior of engineering structures made of the material can be satisfactorily retained. Such idealization can be realized over a representative sample of the composite material. The selection of the dimensions of a representative volume is imperative. The representative volume must be sufficiently large compared to the scale of the microstructure, so that it contains a sufficient number of individual constituents (see, e.g., Neto et al. [13,14]). In addition, it should be entirely typical of the whole composite structure on average (see, e.g., Figure 5). It is assumed that the representative volume contains all microstructural features, i.e.,:  a set of mechanical parameters describing the mechanical properties of constituents,  geometrical parameters characterizing the dimensions of a reinforcement and its location in a representative cell,  a position of a representative cell in a global coordinate system; it may be expressed by Euler's angles in three dimensions or by an angle θ in two dimensions.
For structural scales larger than the representative volume element, continuum mechanics can be used to reproduce the properties of the material as a whole for structural analysis and design without considering the microstructure of the material. Thus, in the macroscopic scale, the microscopically inhomogeneous material can be treated as a homogeneous one, i.e., it is assumed that at each point of the structural body its effective, averaged properties are identical. For a representative volumetric element subject to an imposed macroscopically homogeneous stress or displacement field and no body forces, the average stress and strain components are defined as ( [15][16][17][18]): For structural scales larger than the representative volume element, continuum mechanics can be used to reproduce the properties of the material as a whole for structural analysis and design without considering the microstructure of the material. Thus, in the macroscopic scale, the microscopically inhomogeneous material can be treated as a homogeneous one, i.e., it is assumed that at each point of the structural body its effective, averaged properties are identical.
For a representative volumetric element subject to an imposed macroscopically homogeneous stress or displacement field and no body forces, the average stress and strain components are defined as ( [15][16][17][18]): where V denotes a volume of the representative cell, and x loc denotes the vector of coordinates in the local system of the cell. The effective properties are defined in terms of the relations between the average stresses and average strains over a representative volume [15][16][17][18]: In most structural applications, the structural elements are simplified material models by reducing the general anisotropic stiffness matrix to a simpler form, e.g., to the model of the orthotropic body (Equation (4)). However, the exploitation of the orthotropic model directly in relation (7) isand should be treated only as a hypothesis. The effective properties produced by employing Equation (7) and any simplified model of an anisotropic body always require an experimental verification of results, since the selected assumptions carried out on the stage of Equation (7) may not allow the consideration of certain material characteristics.
For an assumed material model of a two-phase composite, it is possible to evaluate the effective properties with the use of the inverse method by the determination of the concentration matrices in different phases. In the Suquet [15] method, the average strains in each phase are uniquely dependent on the average strains in the representative volume cell: ε where: A f ij U x loc , A m ij U x loc are referred to as concentration matrices for the fibers and the matrix, respectively; and U is a displacement vector in a local (cell) coordinate system. By assuming that the following physical relationship is valid in the local coordinate system: and using the definition of concentration matrices (8a), the following relation holds: Comparing the above equation with Equation (7), the following expression can be obtained: which describes the overall material stiffness matrix. The effective material stiffness matrix of the two-phase composite material is symmetric. A macroscopically homogenized composite material can be introduced (defined) at different levels of the analysis, i.e.,: (i) a laminate or (ii) an individual layer (a lamina). Therefore, the following principle is fundamental in optimization problems considered in this area (see Sections 7 and 8): A laminate (a lamina) is treated as a material continuum, i.e., there are no discontinuities between material phases or individual layers, air bubbles in a polymeric matrix, micro-cracks, etc.
We will find this formalism (commonly called macro-mechanics) useful later in the analysis, when dealing with optimization problems. However, the material modeling of composites refers directly to the total number of independent engineering constants describing mechanical properties. Note also that the number of independent engineering constants is connected with the correct description of structural behavior. On the other hand, it may complicate an optimization problem, especially as the number of constants increases.
In composite material research and technology, a proper characterization of a representative cell, whether with regard to chemical, physical or mechanical properties, is extremely important because it may be the origin of various errors. The prediction of the effective properties for a unit cell proved, and remains, a great challenge. For instance, Figure 6 demonstrates an example of three-dimensional composites where a careful examination of the symmetries of the problem is required.

increases.
In composite material research and technology, a proper characterization of a representative cell, whether with regard to chemical, physical or mechanical properties, is extremely important because it may be the origin of various errors. The prediction of the effective properties for a unit cell proved, and remains, a great challenge. For instance, Figure 6 demonstrates an example of three-dimensional composites where a careful examination of the symmetries of the problem is required. First of all, it is necessary to decide whether to formulate the homogenization problem as a two-dimensional or three-dimensional problem. Figure 7 presents a typical two-dimensional woven composite, whose thickness is small compared with the other dimensions of a unit cell. However, the experiments [19,20] clearly show that variations of weft and warp yarns can be successfully expressed with the use of one coordinate only, but the stiffness matrix has to take at least a two-dimensional form. First of all, it is necessary to decide whether to formulate the homogenization problem as a two-dimensional or three-dimensional problem. Figure 7 presents a typical two-dimensional woven composite, whose thickness is small compared with the other dimensions of a unit cell. However, the experiments [19,20] clearly show that variations of weft and warp yarns can be successfully expressed with the use of one coordinate only, but the stiffness matrix has to take at least a two-dimensional form. increases.
In composite material research and technology, a proper characterization of a representative cell, whether with regard to chemical, physical or mechanical properties, is extremely important because it may be the origin of various errors. The prediction of the effective properties for a unit cell proved, and remains, a great challenge. For instance, Figure 6 demonstrates an example of three-dimensional composites where a careful examination of the symmetries of the problem is required. First of all, it is necessary to decide whether to formulate the homogenization problem as a two-dimensional or three-dimensional problem. Figure 7 presents a typical two-dimensional woven composite, whose thickness is small compared with the other dimensions of a unit cell. However, the experiments [19,20] clearly show that variations of weft and warp yarns can be successfully expressed with the use of one coordinate only, but the stiffness matrix has to take at least a two-dimensional form. The selection of an appropriate homogenization model for heterogeneous composites is based primarily not only on a large number of simplified assumptions dealing with the definition of a unit cell (as described previously), but also on types of prescribed kinematical and static hypotheses characterizing structural deformations. This means simply that the material macro-mechanical model cannot be separated from the structural behavior. They must be compatible in the sense of an accurate formulation of structural deformations. Furthermore, mechanical modeling structure elements composed of fiber-reinforced composites cover both homogenization models and structural models, particularly for composite beams, plates and shells.
To present some of the basic concepts and ideas used in optimization problems, we chose to divide all models (presented in the previous sentence) into four categories, i.e.,: • a plane, orthotropic, two-dimensional (2D) category, • a plane, orthotropic, two-dimensional (2D) category, taking into account transverse shear effects-the stiffness matrix in the form given by Equation (4), eliminating the third row and column in the matrix, • combined two-and three-dimensional approaches, • a three-dimensional (3D) approach used in local formulations of plate and shell theories (layerwise theories).
In the structural analysis of constructions composed of fiber-reinforced plastics, we usually employ the approaches numbered as 1 and 2, and those problems will be thoroughly discussed in the text.It is also worth emphasizing that in formulations 2-4 a set of additional assumptions is required-they deal mainly with the form of transverse shear stress distributions (so-called statical hypotheses) For composite materials reinforced with particles and fibers, a broad review of used material models is discussed in monographs [15,[17][18][19][20][21][22].

Three-Dimensional (3D) Models
In the numerical-based 3D material models, effective properties are evaluated with the use of finite element methods (FEM). They take into account features of the microstructure, such as the fiber architecture. An example of such a model is shown in Figure 8 (Muc [16]). 1/8 of the cylinder (the fibrous reinforcement) embedded in the polymeric matrix constitutes the representative cell. To the representative volume, the prescribed homogeneous displacement-based boundary conditions are imposed (iso-strain method). Usually, the finite element mesh is composed of one to two thousand solid elements. Hence, it is possible to evaluate stress and strain distributions and study the mechanical properties of the assembled composite structure.
 a plane, orthotropic, two-dimensional (2D) category,  a plane, orthotropic, two-dimensional (2D) category, taking into account transverse shear effects-the stiffness matrix in the form given by Equation (4), eliminating the third row and column in the matrix,  combined two-and three-dimensional approaches,  a three-dimensional (3D) approach used in local formulations of plate and shell theories (layerwise theories).
In the structural analysis of constructions composed of fiber-reinforced plastics, we usually employ the approaches numbered as 1 and 2, and those problems will be thoroughly discussed in the text.It is also worth emphasizing that in formulations 2-4 a set of additional assumptions is required-they deal mainly with the form of transverse shear stress distributions (so-called statical hypotheses) For composite materials reinforced with particles and fibers, a broad review of used material models is discussed in monographs [15,[17][18][19][20][21][22].

Three-Dimensional (3D) Models
In the numerical-based 3D material models, effective properties are evaluated with the use of finite element methods (FEM). They take into account features of the microstructure, such as the fiber architecture. An example of such a model is shown in Figure 8 (Muc [16]). 1/8 of the cylinder (the fibrous reinforcement) embedded in the polymeric matrix constitutes the representative cell. To the representative volume, the prescribed homogeneous displacement-based boundary conditions are imposed (iso-strain method). Usually, the finite element mesh is composed of one to two thousand solid elements. Hence, it is possible to evaluate stress and strain distributions and study the mechanical properties of the assembled composite structure.
Various aspects of the application of finite element methods in the homogenization theory are presented, e.g., by Bourget [23] and Guedes, Kikuchi [24].  Various aspects of the application of finite element methods in the homogenization theory are presented, e.g., by Bourget [23] and Guedes, Kikuchi [24].

Two-Dimensional (2D) Models
The macro-mechanical modeling and analysis of a lamina is based on the plane stress assumption. For plies, the elastic behavior is macroscopically quasi-homogeneous and orthotropic, with four independent material parameters: E 1 , E 2 , G 12 , ν 12 . The elastic behavior of a lamina depends on the reference system. The above-mentioned material constants and strain-stress relations are introduced in the material (local) orthogonal coordinate system where the reference axes (x 1 , x 2 ) are parallel and transverse to the fiber direction ( Figure 9)-the so-called on-axis case.
The macro-mechanical modeling and analysis of a lamina is based on the plane stress assumption. For plies, the elastic behavior is macroscopically quasi-homogeneous and orthotropic, with four independent material parameters: E1, E2, G12, ν12. The elastic behavior of a lamina depends on the reference system. The above-mentioned material constants and strain-stress relations are introduced in the material (local) orthogonal coordinate system where the reference axes (x1, x2) are parallel and transverse to the fiber direction ( Figure 9)-the so-called on-axis case. Figure 9. The simplified assumptions used in the analysis of in-plane stress problems for a lamina in the local (material) coordinate system (for laminates, the global coordinate system is rotated around the axis x3).
Using the in-plane stress assumption, it has to be borne in mind that some serious inaccuracies in the mechanical response of a lamina can occur, since the transverse shear stress components are equated to zero, i.e.: σ33= σ13 = σ23 = 0. (9) For plane stress problems, the transverse shear normal strain component 33(3) can be derived from the condition σ33 = 0, namely:      E Q (12) x 3 x 1 x 1 x 2 Figure 9. The simplified assumptions used in the analysis of in-plane stress problems for a lamina in the local (material) coordinate system (for laminates, the global coordinate system is rotated around the axis x 3 ).
Using the in-plane stress assumption, it has to be borne in mind that some serious inaccuracies in the mechanical response of a lamina can occur, since the transverse shear stress components are equated to zero, i.e.: For plane stress problems, the transverse shear normal strain component ε 33 (ε 3 ) can be derived from the condition σ 33 = 0, namely: and finally the stiffness matrix [C] (4) is reduced to the simplified form [5,8]: (11) Q ij (i,j = 1, 2, 6) are stiffnesses reduced to the plane of an individual lamina and are numbered by superscripts in order to distinguish them in a laminate-for instance, the symbol Q (k) ij refers to the k-th layer. They are expressed by engineering constants in the following way: The above stress-stress relations of a lamina are modified in order to take into account thermal, hygrothermal or transverse shear effects. Now, we discuss problems connected with transverse shear deformations of a lamina, since hygrothermal problems will be presented later. In the transverse direction, perpendicular to the material plane (x 1 , x 2 ), a lamina stiffness (strength) is very low, and out-of-plane loading is carried out by a very weak polymeric matrix. To include those effects in the analysis, the stress-strain relations in the k-th layer have to be described with the use of a larger number of engineering constants than previously. This may be done by introducing two additional Kirchhoff's moduli, G 13 and G 23 . Finally, the stresses in the k-th layer are expressed as follows: For thick laminates, as the thickening of a laminate cannot be neglected (e.g.,hygrothermal effects), the transverse shear normal stress component σ 33 cannot be identically equal to zero, and it yields that the full (6 × 6) stiffness matrix given by Equation (4) has to be used in the analysis.

Effective Stiffnesses of Laminates in the 2D Approach
Having formulated the constitutive relations for a lamina composed of a generally orthotropic material, the next step is to derive the constitutive relations accruing from bonding several laminae together. Consider a laminate (hybrid) made of N individual plies. In the global reference system, each layer of a laminate can be identified by its material, its location in the laminate (the integral number k, k = 1, 2, ... N), its thickness t k = z k -z k−1 , and its fiber orientation θ k (Figure 10).

Effective Stiffnesses of Laminates in the 2D Approach
Having formulated the constitutive relations for a lamina composed of a generally orthotropic material, the next step is to derive the constitutive relations accruing from bonding several laminae together. Consider a laminate (hybrid) made ofN individual plies. In the global reference system, each layer of a laminate can be identified by its material, its location in the laminate (the integral number k, k = 1, 2, ... N), its thickness tk = zkzk−1, and its fiber orientation θk ( Figure 10).
As may be seen in the relations (17)- (20), there are symbols indexed by the number of an individual layer (k). They are obtained from Equations (15) by adding the appropriate index to the lamination parameters U 1 , ..., U 7 (it is assumed that the material properties of each layer may be different-a hybrid structure), as well as to a fiber orientation θ in a ply. If the material properties of each lamina are identical, then the relations (17)- (20) can be rewritten in the uniform simplified form. Let us introduce the following notation: (1, cos 2θ k , cos 4θ k , sin 2θ k , sin 4θ k ) z k − z k−1 , The above relation presents 20 quantities in the compact form.Each of them can be obtained by the multiplication of the terms in round and curly brackets. For instance, the quantity V

{B}
(3) is the result of the multiplication of the third term in the round bracket (cos4θ k ) by the second term in the curly bracket ( z 2 k − z 2 k−1 /2) and the summing of all the layers in the laminate. Finally, an arbitrary non-zero component of stiffness matrices (17)- (20) can be expressed with the use of Table 1 and Equation (21).
For plane stress problems, variations of the strain components (ε 1 , ε 2 , ε 6 ) with respect to the z coordinate are prescribed in advance. However, the distributions of the ε 4 and ε 5 strains with the z variable cannot be assumed in an arbitrary way, since they have to satisfy the 3D equilibrium equations for stress components. In the 2D approach it is possible to satisfy the 3D equilibrium equations in an approximate manner only-a so-called statical hypothesis. This may be done by assuming a specific form of constitutive relations for the resultant transverse shear forces and the resultant transverse strain components.For the first-order transverse shear deformation theory (FSDT) in a global reference system, they take the following form: where k 1 and k 2 are transverse shear coefficients. The components of the stiffness matrices can take different values depending on the values of the transverse shear coefficients, and those definitions are different.Vinson [8] proposed to assume the continuous variation of the strain components ε 4 and ε 5 through the laminate thickness in the form of a parabolic function f (z). In this way the stiffness matrix components for the FSDT can be derived in a similar manner as for the 2D plane problem (Equations (17)-(20)), i.e., by the integration over the total laminate thickness: As may be noticed, the two-dimensional macro-mechanical description of laminates is based on two assumptions: • in a lamina, the properties of the fibers and the matrix are smearedout in an equivalent homogeneous material with orthotropic behavior, • the distributions of strains over the total laminate thickness are prescribed in advance, and they have to vary continuously through the total thickness.
In fact, the first one is directly connected with a homogenization theory, whereas the second one constitutes the other form of hypotheses, commonly called kinematical hypotheses. They will be discussed in the next section. For laminates, however, both assumptions have influence on the form of stiffness matrix components. 2D formulations are not sufficient in many problems dealing with composite structural behavior and optimization problems (e.g., delaminations). For such problems, additional assumptions are required to build an appropriate 2D formalism or to use 3D physical relations. Now for an analysis of stress/strain/deformation distributions, the latter approach is typical in the available commercial FE packages.

Strain-Displacement Plate/Shell 2D Relations
In the mechanics of deformable bodies a shell is treated as a three-dimensional body whose metric and shape is identified with a metric and shape of a two-dimensional surface called a reference surface (a mid-surface). The global orthogonal curvilinear coordinate system (the Gaussian system-x, y, z) is placed at the shell's mid-surface, where the z coordinate is normal to the reference surface (Figure 11). From a geometrical point of view, plate theory is a borderline case among shell equations, as the Gaussian curvatures of the shell's mid-surface are identically equal to zero. If a kinematical model of shell deformations is based on the first order transverse shear shell deformation theory (FSDT), a three-dimensional displacement vector is represented in the following form: U 1 (x, y, z) = u(x, y) + zψ 1 (x, y), where: u, v and w are the displacement components of the shell's mid-surface, and ψ 1 and ψ 2 are the rotations of the transverse normals referred to the plane (z = 0)-five independent unknown functions. If the straight lines normal to the mid-surface remain straight and normal to that surface after the deformations (the Love-Kirchhoff hypothesis), then the transverse shear strains ε 4 and ε 5 become negligible, which implies that the rotations of the transverse normals can be expressed in the following way: We hereby obtain the classical Love-Kirchhoff theory (L-K). Thus, in the classical L-K shell theory, the deformations of the shell mid-surface are characterized by three independent parameters (unknown functions). R 1 and R 2 denote the radii of the curvatures at the x and y directions, respectively, and A 1 and A 2 are the Lame parameters, Using Equations (25) or (26), the total strains at arbitrary distance z of the mid-surface are given by: where the components are membrane (mid-surface) strains, and κ are the curvature quantities.

Equations of Motion
For 3D composite structures, let us consider the following functional: where: K denotes the total kinetic energy, Π is the total potential energy and τ means the physical time. The total kinetic energy can be written in the following form: and the total potential energy is defined as follows: Where: W is the potential of internal forces (the strain energy density function): Ω defines the volume occupied by the deformable body. T i is the component of the surface forces in the global reference system. S T is the portion of the surface on which surface forces T are specified. For laminated beams, plates or shells, it is assumed that external loadings are applied on the mid-surface of the structures. U i (x, y, z) denotes the components of the 3D displacement vector defined in a global orthogonal reference system. ρ(z) is the density of a multilayered laminate and can vary with the thickness (z) coordinate: 1( . . . ) is Heaviside's function.
According to the classical macro-mechanical approach for composite structures (the homogenization of ply properties), each component of the functional in Equation (31) is evaluated as the sum of contributions from individual layers. Thus, the laminate strain energy can be expressed in the following way: For beams, plates or shells, the volume element dΩ is represented as: A 1 A 2 dxdydz. For structures considered in the work, the Lame's parameters take the following form: for plates A 1 = A 2 = 1, and for circular cylindrical shells A 1 = 1 and A 2 = R. Using the assumed distributions of displacements U i (x, y, z) (e.g., in the form of Equations (25)) and densities (32) along the thickness direction z, employing physical relations (1) and the general form of strain-displacements relations (24), it is possible to integrate the functional (28) with respect to the z variable. Finally, the functional H is a function of the functions u(x,y), v(x,y), w(x,y), ψ 1 (x,y) and ψ 2 (x,y) and of stiffnesses (17)- (20), characterizing laminate configurations. The derivation of the fundamental set of differential equations can be carried out with the use of the symbolic package Mathematica-one command only.

Static Equilibrium Equations
Finally, for static problems settingthe kinetic energy in Equation (28), we obtain the functional dependent on: the kinematical parameters (functions) characterizing shell deformations at the mid-surface (u, v, w, ψ 1 , ψ 2 ), the Lame parameters A 1 and A 2 and external loadings. The examples of those mathematical operations are demonstrated by Muc [26] both for geometrically linear and nonlinear strain-displacement relations.

Multi-Layered Plates
By setting both Lame parameters to one and carrying out the variations of the functional (30) with respect to the kinematical variables (u, v, w, ψ 1 , ψ 2 ), the equilibrium equations (FSDT) of a laminated multilayered orthotropic plate subjected to an external uniform normal pressure p can be found.
where the stiffness coefficients A 44 , A 45 , A 55 can be derived both from Equations (22) or (23).
A solution of equilibrium equations (34) can be found in an analytical way, as for isotropic plates, using the classical Navier method or the Airy stress functions. For more complicated boundary or loading conditions, solutions can be obtained with the help of the commercial FE packages.

Multi-Layered Circular Cylindrical Shells
Using the first-order shear deformation theory in the global approach and the Euler-Lagrange variational principle, we derive the system of equilibrium equations describing the linear deformations of a circular cylindrical shell. It is assumed that the shell is subjected to a radial pressure p(x,ζ), where x and ζ denote the curvilinear orthogonal coordinates of the shell's mid-surface. After the rearrangements, one can obtain the system of equations in terms of displacements and rotations (u, v, w, ψ 1 , ψ 2 ): where the linear differential operators are defined in the following way: Note that not all operators are given in the explicit form in Equations (36), since some of them are symmetric with respect to the subscripts, i.e., L ij = L ji . The further simplifications of the equations are possible for a specified laminate configuration (both for plates and shells). Now, let us consider a particular form of Equations (36) in the case of: (i) axialsymmetric deformations, (ii) an external loading in the form of axial forces N 0 and a radial normal pressure (they are independent on the longitudinal coordinate x), (iii) a laminate configuration for which [B] ≡ 0 and the stiffness matrix components ( ) 16 = ( ) 26 ≡ 0. The last assumption results in the elimination of the term A 45 ≡ 0. Finally, the equilibrium equations can be expressed in the following form: where: (39) By inserting A 55 = ∞ in Equation (39), the equilibrium equations (38) can be easily transformed to the system of equations valid for the classical Love-Kirchhoff theory. It is worth noticing that the form of equilibrium equations, both for the FSDT theory and the L-K theory, is identical to that for isotropic axially-symmetrical shells (Flügge [27] for the L-K theory). However, for laminated shells, the coefficients in Equations (38) are dependent on a laminate configuration, which may significantly affect the maximum values of buckling or failure (maximum stresses or strains) loads and change a maximal carrying capacity of laminated structures.
For an arbitrary laminate configuration, there is no point in searching for solutions of equilibrium equations given in the general form, (35)- (37), with the use of analytical methods, e.g., an expansion in double Fourier series. Using the Fourier series even for the simplest boundary conditions (e.g., simply-supported b.c.), it is impossible to separate terms corresponding to longitudinal x and circumferential θ variables due to the existence in Equations (35)-(37) of non-zero terms V {...} (4,5) . This complicates the classical optimization analysis significantly, and the correct and accurate evaluation of the objective functions can only be carried out with the help of the FE packages.

First-Ply-Failure Criteria of Laminae (FPF)
For composite structural elements, failure criteria are in fact an extension and adaptation of existing concepts in the mechanics of isotropic bodies (the failure procedures for metallic structures are well established), i.e., the appropriate criteria can be introduced in a linear (maximum stress or strain criteria) or a quadratic form. However, they are mainly based on experimental results for various composites (they are purely empirical), so that there are a lot of different criteria. The major difference between isotropic materials and unidirectional fibrous composite materials is the directional dependence of the strength on a macroscopic scale. All failure theories have a limited area of applications due to a variety of composite materials, their failure mechanisms and different behaviors under combined loading (an interaction of stress components). Failure criteria for composites are formulated both in the stress or the strain space. In general, they have to be consistent with rules of mechanics and principles of mathematics.
Laminate failure criteria are applied on a ply-to-ply basis, and the load-carrying capability of the entire composite is predicted by the laminate theories. A laminate may be assumed to have failed when the strength/strain criterion of any one of its laminae is reached (first-ply-failure). Failure criteria have been established in the case of a layer in its local material coordinate system.
Of all failure criteria available, the following three are considered representative and more widely used in engineering practice and design codes: (i) maximum strain theory, (ii) maximum stress theory, (iii) quadratic (interactive) in the stress space.
Failure criteria are always formulated in the material local coordinate system (on-axis). For composite structures, stress and strain components are evaluated in a global reference system (off-axis), and those values have to be transformed to the on-axis system using classical transformation rules for tensors.
The form of the first two is analogous to that of isotropic structures. They assume no stress/strain interaction. Failure occurs when at least one strain/stress component along one of the principal material axes exceeds the corresponding strain/strength parameter in that direction. In the maximum strain theory the strains are limited in the following way: where as for the maximum stress theory the failure criteria can be written as follows: The ultimate strains are evaluated in the following way (or can be measured): where: Xt longitudinal tensile strength along the fiber direction x1, Xc longitudinal compressive strength along the fiber direction x1, Yt transverse tensile strength perpendicular to the fiber direction x2, Yc transverse compressive strength perpendicular to the fiber direction x2, S shear strength in the plane (x1, x2). In both cases, there are six possible failure modes associated with the sign of stress/strain components for the analyzed plane stress problems. However, the symmetry of shear strength absolute values is expected. Note that due to the non-zero values of Poisson's ratio, failure modes, failure envelopes and failure loads are different for the maximum strength and strain failure theories.
Quadratic (interactive) failure criteria can be expressed in a general form in the on-axis system as follows: For composite materials, a failure envelope represents an ellipsoid elongated along one of its axes. Its shape is dependent on the values of material constants denoted as the symbol F with subscripts. The form of the material parameters uniquely defines the form of the quadratic criterion.
As the most well-known and most often practically used criteria, one can mention the following: the Tsai-Wu criterion (interactive tensor polynomial criterion) [28], where: the Hoffman criterion [29], where: the Puck criterion [30], where: the Tsai-Hill criterion [31] (deviatoric or distortion strain energy criterion), where: For the Tsai-Hill criterion, the subscript U is replaced by the symbols t or c depending on the sign of the stresses.
The quadratic failure criteria do not take into account linear terms σ 6 , because if shear stress is reversed the strength remains the same.
The lamina fails if one of the above criteria is violated. However, the criterion does not give information about how much load can be increased if the lamina is safe or how much it can be decreased if the lamina has failed. We can increase the information given by failure criteria if we use a different variable, which is called the quadratic strength ratio RB. It is defined as the ratio of the maximum (ultimate) load that can be applied {σ} max to the value of the load applied {σ} appl : {σ} max = RB{σ} appl (48) When the allowable stress is reached for the k-th lamina in a laminate (first-ply-failure), then RB=1. If RB > 1.The lamina is safe and the applied load can be increased. RB cannot be less than a unity which has no physical reality.
For linear elastic deformations inserting the relation (48) into Equation (43): (RB) 2 (F xx σ 2 1 + 2F xy σ 1 σ 2 + F yy σ 2 2 + F ss σ 2 6 ) + RB(F x σ 1 + F y σ 2 ) − 1 = 0 For isotropic structures, the strength ratio RB is equivalent to a safety factor. The solution of Equation (46) can be easily found (it must be positive): where: G = F xx σ 2 1 + 2F xy σ 1 σ 2 + F yy σ 2 2 + F ss σ 2 6 , H = F x σ 1 + F y σ 2 (51) We will illustrate the meaning of the strength ratios with a simple example. Figure 12 represents the first-ply-failure envelopes in the strain space for different fiber orientations. Note that the strength ratios in the stress space and in the strain space are the same because the failure criterion is identical. Comparing between the different layers ( Figure 12), whichever layer has the lowest strength ratio would be the first layer to fail, as the load is increasing. Let us analyze a tensile test of a composite structure composed of unidirectional layers oriented at 0 • , 45 • and 90 • . The strength (strain) ratio corresponds to the ratio of the length of the arrow measured to the nearest failure envelope (90 • ) to the length of the line segment measured from the beginning of the arrow indicating the origin of the coordinate system. Of course, since the ratio is less than one, the laminate is unsafe and layers oriented at 90 • will fail for the assumed loading condition. It is worth adding that the value of the strength (strain) ratio is dependent on the form of the loading trajectories.
Although a variety of optimization algorithms exists [32,33] even for rectangular plates, the maximization of FPF loads is not a simple task, due to the complexity of fundamental Equations (34) and non-uniqueness of optimal solutions for different failure criteria (40), (41) or (42). Therefore, the analysis is usually reduced to computations based on the use of the maximum strain theory (40) and the analysis of discrete plate configurations [34][35][36]. For multilayered laminated shells, the analysis of FPF is much more complicated than for plates (see, e.g., the examples discussed in Ref [37]).
For plates and shells with holes, the numerical analysis as well as the experimental data clearly demonstrate that the FPF criteria for notched specimens can be characterized with the use of various functions (see, e.g., Shah et al. [38]). It seems to be much more reasonable to introduce different criteria than classical FPF functions. The objective of the optimization problem can be formulated as the MinMax problem, i.e., to minimize the maximum strain energy density function W (31) along the hole border ∂Ω h . The solutions of such optimization problems are discussed in Refs [39][40][41] for plates and in Ref [42] for cylindrical shells. Although a variety of optimization algorithms exists [32,33] even for rectangular plates,the maximization of FPF loads is not a simple task, due to the complexity of fundamental Equations (34) and non-uniqueness of optimal solutions for different failure criteria (40), (41) or (42). Therefore, the analysis is usually reduced to computations based on the use of the maximum strain theory (40) and the analysis of discrete plate configurations [34][35][36]. For multilayered laminated shells, the analysis of FPF is much more complicated than for plates (see, e.g., the examples discussed in Ref [37]).
For plates and shells with holes, the numerical analysis as well as the experimental data clearly demonstrate that the FPF criteria for notched specimens can be characterized with the use of various functions (see, e.g., Shah et al. [38]). It seems to be much more reasonable to introduce different criteria than classical FPF functions. The objective of the optimization problem can be formulated as the MinMax problem, i.e., to minimize the maximum strain energy density function W (31) along the hole border h   . The solutions of such optimization problems are discussed in Refs [39][40][41] for plates and in Ref [42] for cylindrical shells.

Analysis of Macroscopic (Global) Failure Modes of Composite Structures
The first-ply-failure is one of possible failure global criteria of composite structures. Now, we intend to present criteriaother than FPF that are commonly discussed in the literature.

Concluding Remarks
In the present paper,the attention is mainly focused on the description of homogenization problems, the derivation of the fundamental relations with the aid of the Euler-Lagrange method and the formulation of the first-ply-failure criteria in order to em-

Analysis of Macroscopic (Global) Failure Modes of Composite Structures
The first-ply-failure is one of possible failure global criteria of composite structures. Now, we intend to present criteriaother than FPF that are commonly discussed in the literature.

Concluding Remarks
In the present paper, the attention is mainly focused on the description of homogenization problems, the derivation of the fundamental relations with the aid of the Euler-Lagrange method and the formulation of the first-ply-failure criteria in order to emphasize the influence of the above-mentioned problems on the optimal design of composite materials and structures. In general, my aim was not to solve particular problems but to point out the variety of problems encountered in the design of engineering composite structures.