Plane strain polar elasticity of fibre-reinforced functionally graded materials and structures

This study investigates the flexural response of a linearly elastic rectangular strip reinforced in a functionally graded manner by a single family of straight fibres resistant in bending. Fibre bending resistance is associated with the thickness of fibres which, in turn, is considered measurable through use of some intrinsic material length parameter involved in the definition of a corresponding elastic modulus. Solution of the relevant set of governing differential equations is achieved computationally, with the use of a well-established semi-analytical mathematical method. A connection of this solution with its homogeneous fibre-reinforced material counterpart enables the corresponding homogeneous fibrous composite to be regarded as a source of a set of equivalent functionally graded structures, each one of which is formed through inhomogeneous redistribution of the same volume of fibres within the same matrix material. A subsequent stress and couple-stress analysis provides details of the manner in which the flexural response of the polar structural component of interest is affected by certain types of inhomogeneous fibre distribution.


Introduction
Fibrous composites with either homogeneously or inhomogeneously distributed stiff fibres are increasingly attracting attention and interest, particularly after carbon nanotube fibres were found suitable for inclusion in their constituent phases (e.g., [1]). Despite their low density and nano-meter thickness, carbon nanotubes are known to exhibit remarkably high strength and stiffness, as well as similarly high bending resistance.
Fibre bending stiffness of such a kind of stiff fibres is thus naturally required to be accounted for in modelling and studying the behaviour of relevant composites, regardless of whether fibre reinforcement is distributed in a homogeneous or is some inhomogeneous manner. This requirement becomes particularly important in cases of high fibre concentration (either global or local), where fibre bending resistance gives rise to a couple-stress field. The latter makes the stress field non-symmetric, and endows the composite characteristics of a polar material. It is recalled in this context that the conventional theory of fibre-reinforced materials is built on the simplifying assumption of perfectly flexible fibres [2][3][4], namely fibres that exhibit no bending resistance, and is therefore inherently a non-polar elasticity theory.
The study of polar material behaviour is naturally associated with modelling features falling into the Cosserat theoretical framework [5]. Linearly elastic behaviour of polar fibrous composites may accordingly be attempted through use of either the polar linear elasticity theory proposed by Mindlin and Tiersten [6] for generally anisotropic materials or the linearised version of the theory proposed by Spencer and Soldatos [7] for specific types of polar fibrous composites (see also [8]). It is recalled in this connection that the type of appropriate material anisotropy that fits a relevant boundary value problem is dictated by the specific direction(s) that fibres are aligned to in a fibrous composite.
However, as is also pointed out in [9], there exists no evidence suggesting that the anisotropic version of the Mindlin and Tiersten theory [6] was motivated by potential applications referring to linearly elastic composites having embedded fibres resistant in bending. As a matter of fact, most of the polar linear elasticity analysis detailed in [6] is devoted to the isotropic version of that theory.
Motivated by these observations, the analysis presented in [9]: (i) underlined the principal equations of the Cosserat polar material framework (see also [10]) that provide common ground for the theorie333s proposed in [6] and [7]; (ii) noted the manner in which the linear constitutive equation employed in [6] was obtained through a suitable truncation of the energy expansion proposed by Toupin [11]; (iii) enlarged and enriched the theoretical background through which both theories [6,7] are valid and operate; and, within that enlarged background, (iv) identified their similarities and potential differences without having the intention to either bridge or eliminate the latter.
The principal relevant similarity recorded in [9] refers to the fact that the governing equations of either theory are generally non-elliptic. As a result, the solution to any wellposed boundary value problem, attempted through use of either theory, may be not unique. There are basic historical reasons (see [9]) that prevented Mindlin and Tiersten from noticing this fact in [6], where it is stated that such a potential solution, described by continuous displacements possessing continuous derivatives of all orders, is the unique solution of the implied boundary value problem. . However, Reference [9] has shown that such a solution, which will be termed "the continuous solution" in what follows, is in fact the only possible solution described by continuous displacements possessing continuous displacements of all orders. Due to the observed "lack of ellipticity" of polar elasticity equations, the implied continuous solution may be accompanied by a number of "weak discontinuity" solutions of the same boundary value problem and may thus be not unique. These are solutions described by continuous displacements that possess discontinuous derivatives, and may thus represent micro-scale (fibre-thickness) material failure modes (e.g., [8,12,13]). Such kind of possible solutions are not observable in corresponding problems underpinned by non-polar linear elasticity principles, which always lead to elliptic governing differential equations.
The outlined observations raise immediately a question of whether the prevailing solution of a polar elasticity boundary value problem is the continuous one or some of its possible weak discontinuity counterparts. The task of seeking for an answer to this question is of paramount practical importance in structural analysis applications.
Such a challenging task may well depend on the particular polar elasticity problem of interest. Moreover, it seemingly requires some analytical and/or numerical/computational comparison of all relevant weak discontinuity solutions among themselves, as well as against their common continuous counterpart. In fact an appropriate comparison may also be required of the stored energy levels reached by all possible solutions involved. The need becomes thus evident for the derivation of relevant continuous and/or weak discontinuity solutions to a number of relatively simple or more difficult boundary value problems, with the aim to reach afterwards a stage that makes the implied comparisons possible.
The present study is considered as an immediate continuation of an initial, relatively simple step made already in that direction [14], in the sense that it complements the latter in the search for continuous solutions to the plane strain bending problem of a simply supported, linearly elastic rectangular strip reinforced by a single family of straight fibres resistant in bending. While [14] dealt with the case of either homogeneous or layer-wise inhomogeneous (laminated) strips, this communication considers the more general case of material inhomogeneity due to continuous through-thickness variation of the fibre-reinforcement.
The elastic strip of interest may be regarded as a rectangular beam made of functionally graded material (FGM) having unit width, or as the cross-section of a corresponding rectangular plate having infinite extent in the out of plane direction. The latter representation provides direct connection with the relevant, non-polar elasticity problem considered and solved by Pagano [15] but, here, the implied bending stiffness of functionally graded fibres furnishes the strip with polar material properties. Moreover, material inhomogeneity features in the analysis through the variable form attained by the coefficients of the corresponding set of Navier-type partial differential equations.
With the help of Appendix A, Section 2 thus provides a proper mathematical description of the plane strain state of polar, linearly elastic structures reinforced in a functionally graded manner by a single family of straight fibres resistant in bending. For simplicity, this description is based on the restricted version of the polar elasticity theory presented by Spencer and Soldatos [7]. This version of the theory (see also [8]) involves only a single elasticity modulus of fibre bending resistance and, as soon as certain additional conditions are met [9], can establish connection with the theory of Mindlin and Tiersten [6].
Section 3 formulates the aforementioned bending problem of a simply supported prismatic beam (or rectangular plate cross-section). Moreover, with use of a second Appendix, it employs a suitable semi-analytical mathematical method (e.g., [16][17][18][19][20][21]), provides information that underpins its computational efficiency, and finalises the solution of the corresponding Navier-type equations. With the help of Appendices C and D, Section 4 connects afterwards the present problem of interest with its homogeneous polar elasticity counterpart [14]. This connection enables a homogeneous fibre-reinforced component [14] to be regarded as the source of a set of equivalent functionally graded structures, each one of which is made through inhomogeneous redistribution of the same volume of fibres within the same matrix material.
Three different types of such inhomogeneous fibre redistribution are thus selected in Section 4, and are employed afterwards in Section 5, in the discussion of the numerical results presented there. That Section thus examines in detail the manner in which each of the employed types of inhomogeneous fibre distribution affects the flexural response of the composite structure. Finally, Section 6 summarises the main conclusions drawn and outlines directions on the manner in which the search for identification of corresponding weak discontinuity solutions should be contacted.

Theoretical formulation
A linearly elastic fibre-reinforced plate has finite length, L, in the x direction, infinite extent in the y direction, and finite thickness, h, in the z direction of a Cartesian co-ordinate system The plate material has embedded a single family of fibres which are parallel to the x-axis, can resist bending, and are distributed in the z-axis direction in a continuous, functionally graded manner. The plate is subjected to external loading that justifies plane strain response, in the sense that the displacement component in the y direction is zero while the other two displacement components, as well as all remaining physical quantities, are independent of the co-ordinate parameter y. In the usual manner, the plate cross-section can thus be considered as a two-dimensional elastic strip or as a prismatic beam having length L, thickness h and unit width ( Figure 1). In this context, relevant terminology of prismatic beams is also employed in what follows.
The through thickness inhomogeneous distribution of the fibres is regulated by controlling their volume fraction, ( ), which requires from the material properties of the structural component to be known functions of the z co-ordinate parameter. Every material property, P(z) say, of such a functionally graded fibrous composite is usually expected to obey the mixture law where V m (z) is the volume fraction of the matrix phase, while P f and P m represent the corresponding constant material property of the fibre and the matrix phase, respectively. It is pointed out that the inequality conditions noted in (1) are imposed on the ground of evident theoretical arguments that hold true regardless of the particular form of ( ) or, equivalently, ( ). In this context, the denoted upper limit of the fibre volume fraction, namely ( ) = 1, is in principle possible only in cases that fibres are assumed perfectly flexible and, having no thickness, can therefore fill in completely the entire volume of the composite. However, fibres do have thickness in practice and, due to the structural architecture of the fibrous composite, leave among them gaps which are filled in with matrix material even in parts of the composite that fibres are distributed very densely.
A more realistic approach thus requires introduction of a maximum fibre volume parameter, say, such that This additional condition does not need to be discussed further at these early stages of the problem formulation. However, it is reconsidered and discussed later, in Sections 4 and 5, where determination of a value for becomes part of some specific applications. In the implied plane strain state, the average fibre and matrix concentrations of the composite are defined as follows: The particular case of a homogeneous fibrous composite, where the fibre volume fraction is constant, is thus characterised by the relationship ( ) = < > for all z. If the fibres resist bending and < > is adequately high, say 40% to 60%, then the fibre response to mechanical loading generates considerable couple-stress and non-symmetric stress (e.g., [7][8][9]14]). In the case of FG fibrous composites with relatively low < >, creation of a couplestress field is still possible locally, namely in specific parts of the composite where ( ) anticipates high fibre concentration.
The stress and couple-stress components that contribute actively in plane strain equilibrium are shown schematically in Figure 2 (see also [7,14,22]). The symmetric part of the stress tensor is given by the standard form of the generalized Hooke's law, which in the present, plane strain case acquires the form where the appearing linear strain components are = , , 2 = , + , , = , .
Here, U(x,z) and W(x,z) are the displacement components along the axial and transverse coordinate direction, respectively, and a comma denotes partial differentiation with respect to the indicated co-ordinate parameter(s). The elastic moduli appearing in (4) vary in the transverse direction in accordance with the mixture law (1), namely As the matrix phase is naturally considered isotropic, the following relationships are assumed valid: where λ and μ are the constant Lamé moduli of the matrix material. In this context, Appendix A describes an alternative manner in which the elastic moduli of the matrix and the fibre phases can be related, and thus lead to the determination of their Cij counterparts appearing in (4). The antisymmetric part of the stress tensor is defined as follows: where the only non-zero couple-stress component met in this plane strain problem (e.g., [7,14,22]), namely acts in the manner shown in Figure 2, and represents the fibre curvature. Unlike Cij which have dimensions of stress, the fibre bending modulus d f has dimensions of force. Like Cij though, this is also expected to obey the mixture law (1).
However, unlike the fibre phase, the isotropic matrix phase does not contribute to the bending stiffness of the fibrous composite, and, as a result, the second term in the right-hand side of the corresponding expression (1.a) is zero. Hence, in line with previous relevant studies [14,22], where material homogeneity enabled the relevant constant value of d f to be considered as a product of the form 11 , the fibre bending modulus attains here the through thickness variable form where the intrinsic material length parameter l is connected with the fibre thickness. In this manner, l = 0 represents cases of non-polar material behaviour, where fibres are perfectly flexible and the subsequent absence of couple-stress ( , = [ ] = 0) enables the stress tensor to attain its conventional symmetric form (4).
When l is non-zero, the non-zero shear stresses are unequal, so that In the absence of body forces, the equilibrium equations thus acquire the form which, after appropriate use of equations (4)- (9), lead to the Navier-type equations The outlined polar elasticity formulation is general, in the sense that it applies to all cases that a relevant FGM fibrous composite exhibits plane strain behaviour. For analytical purposes, it is found convenient to rearrange equations (13) into the following matrix form:

Cylindrical bending of a simply supported plate
Attention is now confined into the particular case that deformation is due to external application of the lateral boundary tractions The externally applied transverse load, q(x), is considered known and can, therefore, be represented in the following Fourier-type sine-series form: It is further assumed that the longitudinal ends of the plate cross-section or prismatic beam (x = 0, L) are subjected to the following set of homogeneous boundary conditions:  which is consistent with the symmetries of simply supported boundaries. In the particular case of a homogeneous fibrous composite, where V f and V m are both known constants, the problem of present interest thus reduces naturally to its polar elasticity counterpart studied in [14]. The simple support boundary conditions (17) are satisfied exactly by the following choice of a displacement field: where the functions ( ) and ( ) are to be determined. Expressions (18) represent a potential solution to the described boundary value problem when the external loading is identical with a single term of the series expansion (16), namely The linearity of the described boundary value problem, combined with the superposition principle of relevant solutions, makes it then sufficient for someone to look only for a solution of the particular case in which the external load is given according to (19), with m being an arbitrary positive integer. Upon inserting (18) into (14), the latter equation is transformed into a fourth-order set of simultaneous ordinary differential equations (ODEs) with variable coefficients. This can be expressed as follows: where, ( , ) = [ 1 + 2 2 + 3 ( 4 + 11 ) + 5 − 4 2 + 7 6 + 8 + 9 2 + 12 ], Due to the variable form of ( ), the appearing coefficients, namely 1 = −ℎ 2 2 , 2 = 55 / 11, 3 = 55, / 11, 4 = ℎ ( 13 + 55 )/ 11 , 5 = ( ℎ/ )( 11 / 11 ), 6 = 33 / 11 , Solution of (20) is here achieved with the use of a semi-analytical method, which considers that the inhomogeneous polar material strip of interest is essentially made of an infinite number of fictitious layers having infinitesimally small thickness and constant material properties. As computational practice requires use of a finite number of such fictitious layers (see Figure 1), the larger the number of those fictitious layers considered the nearer the obtained numerical results approach their exact elasticity counterparts.
The implied "fictitious layers method" was initially introduced for the solution of nonpolar linear elasticity problems dealing with the dynamic response of isotropic cylindrical components [16]. In such problems, it is the geometry rather than the material inhomogeneity of the structure that spreads variable coefficients into the governing differential equations. The method has since been applied successfully to both static and dynamic studies of homogeneous and laminated composite components of cylindrical geometry (e.g., [16][17][18][19][20] and relevant references therein), and is proven capable to provide asymptotically identical results to those based on potential or existing exact elasticity solutions.
Moreover, the numerical stability and the rate of convergence of the method are found in practice superior to those of corresponding analytical solutions based on power-series methods (e.g., [20]), where computational practice still requires some suitable finite term truncation of ultimately infinite series solutions, and, hence, does not avoid the concept of an approximation. More recently, the applicability of this fictitious layer method has successfully been extended towards solution of relevant structural mechanics problems that involve even doubly curved functionally graded structural components [21].
Description of the solution thus obtained is facilitated by initially converting (20) into the following, equivalent set of four first-order linear ODEs with variable coefficients: The implied solution then continues by resembling its counterparts described in [16][17][18][19][20]. For self-sufficiency of this communication, further details are briefly presented in Appendix B.

Application for selected forms of the fibre volume fraction
As is pointed out in Section 2, even in parts of the composite where fibres are distributed very densely, the fibre structural architecture leaves gaps which are naturally filled in with matrix material. The inequality conditions noted in (1) should accordingly be refined through use of the more realistic inequality (2), provide that proper consideration of the fibre-scale structure can enable determination of the refined upper bound parameter . This may be achieved with use of some appropriate representative volume, or area elements (RVE) of the fibre distribution pattern [23]. Nevertheless, any -value thus obtained depends on the chosen RVE discretisation, and may therefore be not unique.
Appendix C thus demonstrates the manner in which rectangular or triangular RVEs of the kind implied in Figure 3c can be used as reasonably simple examples in the present problem of interest, where the direction of the considered uniaxial family of fibres in normal to the depicted yz-plane (see also Figure 1). The two different values of maximum fibre volume fraction thus determined in (C.5) and (C.6) are here conjoined as follows: max 0.785 for rectangular RVEs, 4 0.907 for triangular RVEs. 23 4.1 Particular case: Homogeneous composites [14] In order to deal with applications of the outlined analysis, connection is initially established with the corresponding study detailed in [14] for corresponding homogeneous fibrous composites. It is recalled in this context that the effective material properties of the homogeneous fibrous composite employed in [14] are as follows: The analysis is evidently capable to consider homogeneous fibre distributions by using appropriate constant values of the fibre volume fraction and, upon taking (2) and (3) into consideration, it thus requires Let us, for instance, consider the choice which refers to a homogeneous fibrous composite whose volume consists 50% homogeneously distributed fibres and 50% matrix material. Upon inserting into (A.3), and making also use of (3), the mixture law (1) reveals that the effective material properties of the corresponding homogeneous composite are those detailed in (26). With use of (A.2), the holding relationships between the elastic moduli of the corresponding fibre and matrix phases are then found to be as follows: It can then readily be verified that, in this particular case that the fibrous composite of interest is homogeneous and possesses effective material properties of the kind described in (26), the present analysis produces identical displacement and stress distributions to those detailed in [14] with its first iteration (N = 1). Further iterations are unnecessary in that case, as they naturally return the same numerical results.
It is emphasised that the outlined verification of the present analysis is still possible for constant choices of that differ from (28), as soon as the values of the constants (29) and, subsequently, of the ratios (30) are modified in a manner that enables the mixture law (1) to yield again to the effective elastic moduli (26). A couple of specific, additional relevant cases are in fact identified in Section 4.3 below, in connection with the form (35) of possible inhomogeneous fibre distribution.

Functionally graded, linear redistribution of the fibres
The connection established with the homogeneous fibrous composite considered in [14] is now exploited by considering the following pair of linearly inhomogeneous fibre distributions: As either of these return they both represent corresponding inhomogeneous composites consisting 50% fibre and 50% matrix material. A schematic representation of the volume fraction of these inhomogeneous fibre distributions is depicted in Figure 4, along with their counterparts that represent the homogenous composite described already in Section 4.1 (Fig. 4a). Both (31) and (32) are evidently fibre distributions which vary linearly through the thickness, and are non-symmetric with respect to the middle plane of the composite beam. It is evident that (31) represents a top-stiff fibrous composite while (32) corresponds to a bottom-stiff such.
The relevant numerical results presented in Section 5 below refer to inhomogeneous composites whose effective material properties are evaluated with use of the mixture law (1), after each of (31) and (32) is inserted into (A.3). This process requires also use of (29) and (30), so that the resulting inhomogeneous composite is thought of as formed by a relevant redistribution into the same matrix of a same volume of fibres (50%) possessing the material properties (30). It can indeed be readily verified that, in both cases, the obtained through thickness average elastic moduli are still in exact agreement with the effective material properties (26) of the homogeneous fibrous composite employed in [14].

Symmetric, piece-wise linear redistribution of the fibres
The last fibre distribution of present interest is associated with a class of inhomogeneous fibrous composites whose volume fraction varies symmetrically with respect to the middle plane of the composite (see Appendix D). This class is described as follows: where is some real positive constant. In this case, the fibre volume fraction increases in a piece-wise linear manner with increasing the distance from the middle plane (see Figure 4d). As (0) = 0 and (ℎ/2) = (−ℎ/2) = /2, (35) refers to a fibrous composite graded in a manner that maximum fibre volume fraction is attained at both lateral planes. As is also noted in Appendix D, (35) yields the average fibre volume fraction (33) if = 2.
However, as is shown next with the help of Appendix C, the types of fibre-scale structure considered there make (33) incompatible with the fibre distribution (35). This is because connection of (35) with the maximum volume fractions noted in (25) gives, respectively, the following maximum value of the positive constant : .
A considerable part of the numerical results presented in the next Section thus refers to the polar mechanical response of this inhomogeneous fibre-reinforced composite which, along with employing the material properties (38), implies further that = 1.814 in (35).

Numerical results and discussion
All numerical results presented and discussed in this Section are obtained by setting m = 1 in (19). These results are presented in a non-dimensional form, through use of the following dimensionless displacement and stress parameters ̅ = < > By virtue of (33), these are seen equivalent to their counterparts employed for the corresponding case of a homogenous fibrous composite in [14], where, however, an evident typographical error is noticed in the couple-stress non-dimensionalisation. The evident symmetries that (19) imposes along the x-direction imply that the magnitude of displacements, stresses, and couple-stress have identical through-thickness distribution at x/L and 1 -x/L. Numerical results are accordingly presented for the left half of the beam only. All rectangular beams considered for the results shown next possess the same span ratio with their homogeneous counterpart studied in [14], namely For a natural connection with [14], the same notation, namely is used for the non-dimensional intrinsic material parameter that refers to fibre thickness. In this regard, a note is made of the fact that this parameter should not be misinterpreted as denoting the Lamé modulus employed in (7) and (A.1).
As is connected with the fibre thickness and, hence, cannot exceed the beam thickness, λ acquires naturally the upper bound value noted in (41) only if = ℎ. However, connection of with the fibre thickness is here refereed to only as an example of the manner in which one can handle the aforementioned dimensions difference between the fibre bending modulus, d f , and the conventional elastic moduli met in non-polar elasticity.
If, for instance, one accepts that fibres are approximately arranged through the beam thickness in the form of representative volume elements described in Appendix C, the estimated upper bound of λ may further be reduced considerably, and/or even be related to the It is pointed out though that, in view of (10), (41) leads essentially to the following reparametrisation the fibre bending stiffness modulus: This relationship shows that, although useful on physical grounds, λ is not necessarily the most influential parameter for a proper determination of . In fact, determination of in a structural component should still be based on potential experimental work and observation, precisely as happens with the determination of conventional elastic moduli.
By setting = 0 in (31) or (32), it is thus made initially sure that the present analysis gives identical numerical results, and is thus in complete agreement with its counterpart presented in [14]. This confirmation enables next consideration and study of corresponding numerical results that refer to inhomogeneous relevant composites having fibre volume fraction of the type (31) and (32) with ≠ 0, or (35) with = 1.814.

Through thickness displacements distributions
For several different values of the fibre inhomogeneity and the fibre bending stiffness parameters, Tables 1 and 2 present the non-dimensional value of the in-plane and the transverse (flexural) displacement, respectively, obtained at selected points through the thickness of a top-stiff beam. To a considerable extent, these results are susceptible to comparison with their counterparts presented for in Tables 1 and 2 of [14], respectively, for corresponding homogeneous fibrous composites.
In line with the relevant trend noted in [14], Table 2 thus confirms that the magnitude of the flexural displacement decreases with increasing the fibre bending stiffness parameter, λ, due to the additional flexural stiffness provided by the fibre bending resistance. However, it is seen here further that the magnitude of the flexural displacement decreases further with increasing the inhomogeneity parameter, ε. This is because, by increasing ε, the bending stiffness of the beam is increasing near the top lateral boundary where the external load is applied.
It is recalled on the other hand that the results presented in Tables 1 and 2 of [14] show that, in the case of a homogeneous beam (ε = 0), the in-plane displacement is always at least an order of magnitude smaller than its flexural counterpart. However, upon increasing the non-zero value of ε, the increasing material inhomogeneity affects the existing local coupling between bending and extension to such an extent that the magnitude of U becomes comparable to that of W .
It is then not surprising that the values of U shown in Table 1 differ from those of their counterparts presented in [14] even for ε = 0.01. In fact, for ε = 0.05 the values of U are already comparable with their W counterparts (Table 2). Moreover, for ε = εmax = 0.812, which is the maximum value assigned to ε when fibre scale structure is designed with use triangular RVEs, the magnitude of the in-plane displacement parameter exceeds that of W , at least within the adopted region of the λ-variation. It is pointed out that, as all numerical results shown in this study refer to the left half of the beam, the minus sign associated with almost all numerical values shown in Table 1 implies that the beam deformation creates a predominantly tensile in-plane displacement.
Analogous conclusions may be drawn by observing and comparing the numerical results tabulated in Tables 3 and 4 for corresponding U -and W -values of a bottom-stiff inhomogeneous beam. The magnitude of displacements is again decreasing with increasing the value of the fibre bending stiffness parameter, λ. However, the sign of almost all numerical values shown in Table 3 reveals that the beam deformation creates now a predominantly compressive in-plane displacement.
Strong local inhomogeneity effects, of the type observed previously in Tables 1 and 2, have now emerged mainly at the bottom part of the beam. It is instructive for someone to observe that for ε = 0.05 the magnitude of the in-plane displacement (Table 3) is again comparable with its flexural displacement counterpart (Table 4). Moreover, for ε = εmax = 0.812, the former parameter exceeds the latter to such a substantial degree, that the deformation seems in this case to take mainly place through in-plane extension rather than flexure. Nevertheless, Tables 1-4 suggest that, in general, top-stiff beams suffer smaller flexure and, therefore, may generally be considered stronger than their bottom-stiff counterparts at the same value of the inhomogeneity parameter, ε.
The observed in-plane deformation dominance seems to increase with increasing ε to an extent that affects substantially the detailed features of relevant stress distributions. This observation is briefly discussed later in Section 5.3, which illustrates the influence that the increasing value of ε exerts on the bending stress distribution observed within both top-and bottom-stiff beams.
Under these considerations, the corresponding non-dimensional displacement results shown in Tables 5 and 6 suggest that beams with through thickness symmetric fibre distribution are similarly strong. Indeed, the magnitude of the flexural displacements shown in Table 6 are comparable with their counterparts shown in Table 2 for the largest value of the fibre inhomogeneity parameter, ε = εmax = 0.812, at least within the top half of the beam.
Moreover, while the top-stiff beam has higher average fibre volume fraction, the difference observed between corresponding numerical results shown in Table 2 and 6 is decreasing at the top part of the beam with increasing the fibre bending stiffness parameter

Through thickness couple-stress and shear stress distributions
As transition from non-polar to polar material behaviour is caused by the emerging couplestress field, immediate attention is next directed towards the influence that couple-stress creation exerts on the shear stress components, giving thus rise to non-symmetric stress. Corresponding numerical results that show the manner in which normal stresses are affected are also presented and discussed afterwards, in Section 5.3.
In this context, Figures 5 and 6  For different values of the fibre bending stiffness parameter, , Figures 7 and 8 illustrate next the through-thickness distribution of the shear stress ̅ at the left end of a topand a bottom-stiff beam, respectively. In line with [14], all depicted distributions satisfy the zero shear traction boundary conditions imposed on the upper and lower surface of the beam. Due to the relatively small value of the material inhomogeneity parameter (ε = 0.1), the depicted curves do not diverge substantially from their counterparts shown in Figure 6 of [14]. However, they have all lost their largely symmetric form observed in [14] with respect to the beam middle axis, while their highest magnitude is moved towards the direction of increased fibre reinforcement; namely, upwards for the top-stiff and downwards for the bottom-stiff beam. As the beam becomes stiffer with increasing , that highest magnitude of ̅ decreases and moves naturally towards zero. It comes, however, a little as a surprise that, in the case of the bottom-stiff beam ( Figure 8) and for the relatively large value  = 0.009 of the fibre thickness parameter, the relatively small value of that highest ̅ -magnitude changes sign, along with the sign of the whole ̅ -distribution. Figures 9 and 10 show next the through-thickness distributions of the stresses ̅ that correspond to ̅ -distributions illustrated in Figures 7 and 8 Due to the relatively small value of the fibre inhomogeneity parameter (ε = 0.1) the depicted ̅ -distributions present again similarities with their counterparts shown in Figure 5 of [14]. Nevertheless, in almost all cases, the highest magnitude of the ̅ -value moves again towards the stiffest part of the inhomogeneous structural component. An exception to this trend is again observed in the case of the bottom-stiff beam ( Figure 10) where, for the relatively large value  = 0.09 of the fibre thickness parameter, the ̅ -distribution reverses hollows and attains highest magnitude on the bottom lateral boundary.
Another remarkable observation relates to the fact that, like their counterparts depicted in Figure 6 of [14], all ̅ -distributions shown in Figures 9 and 10 intersect at a certain pair of material points located at the vicinity of z/h = ±0.3. At those points, the value of ̅ thus seems independent of the fibre thickness parameter, , although it evidently still depends on the fibre inhomogeneity parameter, ε. At present, there seems no obvious explanation to this effect, which is apparently due to the manner in which the couple-stress influences the values of ̅ .
Under these considerations, Figures 11 and 12 depict the through thickness ̅ -and ̅ -distributions, respectively, predicted at the left end of a beam reinforced in the symmetric, piece-wise linear manner (35). Remarkably, changes of the fibre bending stiffness parameter, , do not seem to influence notably either of these shear stress distributions. The considerable similarity observed between the ̅ -distributions depicted in Figure 11 and their ̅ counterparts shown in Figure 12 is thus not surprising.
In fact, the principal difference between corresponding results demonstrated in those Figures is that all ̅ -distributions ( Figure 11) attain, naturally, a zero value at the top and bottom lateral plane. Like their counterparts depicted previously in Figures 9 and 10, all different ̅ -curves shown in Figure 12 pass again through a certain pair of material points, which are now moved at the vicinity of z/h = ±0.4. Moreover, the lack of fibre-reinforcement on the middle-axis has apparently made z/h = 0 a third point of intersection for all ̅and ̅ -curves depicted in Figures 11 and 12. It is thus observed that, due to low local fibrereinforcement, the stress state is nearly symmetric within a certain material band that surrounds the beam middle-axis.  [24]. The remaining curves then show the influence that fibre bending stiffness exerts on those known distributions upon gradually increasing the value of . Figure 13 thus makes it clear that fibre bending stiffness has marginal effect on the transverse normal stress distribution.

Through thickness normal stress distributions
However, Figure 14 reveals that, upon increasing , the increasing resistance of the beam lowers the magnitude of ̅ and, hence, decreases the influence that the depicted, wellknown boundary layer effect of the ̅ -distribution exerts on the strength of the structure. In fact, for  as small as 0.009, that influence is decreased to such an extent that the value of the axial normal stress might be felt notable only within a particularly narrow layer near the beam lateral boundary. Still though, every curve shown in Figure 14 evolves about a pivotal point, located at the vicinity of ( ̅ , ) = (0, 0), in a manner that divides the corresponding bending stress distribution into a compressive and a tensile part (top and bottom half of the beam, respectively).
In the light of these observations, Figures 15 and 16 present next evidence of the fact that combined action of fibre bending stiffness and material inhomogeneity (ε = 0.1) has still marginal effect on the ̅ -distribution of a top-and a bottom-stiff beam, respectively. However, the same is not true for the corresponding ̅ -distributions.
A search for the effect that combined action of fibre bending stiffness and material inhomogeneity has on the ̅ -distribution is facilitated by initially considering that fibres are perfectly flexible ( = 0) and varying only the value of the inhomogeneity parameter ε. In this context, Figures 17 and 18 demonstrate the manner in which the "blue" ̅ -distribution associated in Figure 14 with  = 0 evolves with increasing inhomogeneity in a top-and a bottom-stiff beam, respectively. A thorough study of the results presented in Figures 17 and  18 makes afterwards clearer the corresponding results depicted in Figures 19 and 20 for corresponding inhomogeneous beams with embedded fibres resistant in bending ( = 0.06) Figure 17 thus reveals that, upon increasing ε, the observed ̅ -distribution is initially transposed to the right. This is due to the dominance the ε exerts on its linearly dependent elastic moduli for small amounts of inhomogeneity. Hence, upon increasing ε within a narrow interval of relatively small values, the corresponding change of the elastic moduli suffices to "push" the depicted curve to the right, to an extent that soon turns the whole ̅ -distribution completely tensile. Nevertheless, beyond the value ε = 0.083, which is still relatively small, the inhomogeneity difference between the top and the bottom parts of the beam becomes very dominant. Upon increasing ε further, the ̅ -curve thus changes shape and, while still moves to the right and hence stays tensile, reveals that it is the upper, rather than the bottom part of the top-stiff beam that bares most of the externally applied loading.
Eventually, at about ε = 0.25, the top part of the beam becomes so stiff that, while the value of the tensile bending stresses start to decrease at the lower part of the beam, the value of ̅ approaches a maximum on the top lateral plane. Further increase of the ε-value and, hence, of the stiffness of top part of the beam lowers the observed tensile bending stresses throughout the beam thickness. In is instructive in this regard to note that the ̅ -curve drawn in Figure 17 for ε = 0.4 is essentially transposed to the left when the inhomogeneity parameter is increased up to ε = 0.5708, or further up to ε = 0.812. It is recalled that (34) associates 0.5708 or 0.812 with the maximum value that ε may attain when the fibre-scale structure is modelled with rectangular or triangular RVEs, respectively. It is thus observed that, upon increasing ε towards its maximum value, the decreasing tensile value of ̅ observed near the bottom boundary of the top-stiff beam is naturally moving towards zero. Figure 18 demonstrates that the bending stress distribution of a bottom-stiff beam with embedded perfectly flexible fibres responds in an analogous manner. One of the evident principal differences with the trends observed in Figure 17 is the fact that, upon increasing material inhomogeneity, the ̅ -distribution that corresponds to  = 0 in Figure 14 moves towards the left, and thus soon turns completely compressive. Moreover, the change of shape that the ̅ -curve observes for higher values of ε suggests that, naturally, it is now the bottom, rather than the upper part of the bottom-stiff beam that bares most of the loading. Finally, the aforementioned observations, associated in Figure 17 with the top and the bottom plane of a top-stiff beam, are naturally seen in Figure 18 associated with the bottom and the top plane, respectively, of a bottom-stiff beam.
In the light of these observations, Figures 19 and 20 illustrate next the manner in which the "yellow" ̅ -distribution, formed in Figure 14 by setting  = 0.06, evolves with increasing inhomogeneity of a top-and a bottom-stiff beam, respectively. To a considerable extent, these results show substantial quantitative similarity with their counterparts depicted in Figures 17 and 18 for corresponding beams having embedded perfectly flexible fibres. However, and in close agreement with all previously observed trends, the extra bending stiffness added now on the functional graded beam lowers significantly the magnitude of the observed bending stresses. The two final Figures, namely Figures 21 and 22, depict the distribution of the normal stresses ̅ and ̅ , respectively, at the mid-span of a fibre-reinforced beam subjected to the symmetric, piece-wise linear fibre reinforcement (35). In line with corresponding results depicted in Figures 11 and 12 for corresponding shear stress distributions, both stress distributions attain a nearly symmetric form, with respect to the middle axis, within a narrow band of weak local fibre-reinforcement. Regardless of the value of the fibre bending stiffness parameter, that symmetry gradually fades outside that band, which surrounds the middle axis.
The value of the bending stress at the top boundary of the beam, where the external load is applied, thus becomes naturally bigger from its bottom boundary counterpart (see Figure 21). However, like the top-stiff beam (see Figures 17 and 19 for ε ≠ 0), the bending stress is always tensile throughout the beam thickness. Unlike the corresponding results shown in Figures 17 and 19 though, the imposed lack of local fibre-reinforcement at z/h = 0 enables the beam middle axis to remain almost free of bending stress. It is also noticeable in this regard that, regardless of the value of the fibre bending stiffness parameter, the through thickness shape of the ̅ -distribution resembles closely the form (35) of the corresponding fibre volume fraction.
On the other hand, Figure 22 reveals that the weak reinforcement observed around the beam middle axis and, hence, practically the negligible influence of ̅ is locally compensated by a sharp jump of the transverse normal stress, ̅ . The latter is seen positive and, therefore, tensile within the aforementioned narrow band, as well as within the bottom part of the beam. Unlike its top-and bottom-stiff counterparts, which are negative and therefore compressive throughout the beam thickness (see Figures 15 and 16), the distribution of ̅ depicted in Figure 22 is compressive at and near the top beam boundary, where the external load is applied, but turns tensile within the aforementioned band of weak fibre reinforcement. It then remains tensile in the bottom part of the beam, where it decreases and becomes finally zero on the unstressed bottom beam boundary.

Conclusions
This study aims initially to promote the need for extension into the regime of polar material response of fibrous composites of relevant non-polar linear elasticity solutions. Namely, existing solutions of well-posed boundary value problems represented by continuous displacements having continuous derivatives of all orders. As non-polar elasticity of fibrereinforced materials assumes that fibres embedded in a structural component are perfectly flexible, the implied solution extensions will offer substantially better understanding of the behaviour of composites reinforced by stiff fibres, such as carbon nanotubes, that exhibit considerable bending resistance.
In serving the first of these aims, this study continued an initial relevant extension [14] of the well-known Pagano's non-polar plane strain elasticity solution [15], by considering that the implied stiff fibres are redistributed within their matrix in an inhomogeneous, functionally graded manner. Like in [14], the implied solution extension was based on the restricted version of the polar elasticity theory presented by Spencer and Soldatos [7], namely a theory that involves only a single elasticity modulus of fibre bending resistance. That extra elasticity modulus is associated with the global response of the fibrous composite, rather than with the response of individual fibres, but its involvement enables the theory to make use of an intrinsic length parameter that relates to an average fibre thickness.
By setting that intrinsic length parameter and, therefore, the fibre thickness equal to zero, the present theory and analysis reduce naturally to their conventional, non-polar elasticity counterparts. It becomes thus understood that the content and the results of this article are useful even in cases of non-polar material response, where the fibres embedded in a relevant functionally graded fibrous composite are perfectly flexible. In this context, the parametric studies performed in Section 5, along with their counterparts presented in [14], enable better understanding of the influence that fibre bending resistance exerts on the plane strain behaviour of the implied class of fibrous composites, provided that the obtained, continuous solution of the boundary value problem (Section 3) prevails over potential weak discontinuity solutions.
In this connection, it is reemphasised that this communication aims further to make it wider known that, unlike its non-polar linear elasticity counterpart, a corresponding fully continuous polar linear elasticity solution is not necessarily the unique solution of the respective boundary value problem. Due to the lack of ellipticity that linear elasticity equations suffer outside their non-polar material regime, the implied polar elasticity extension of a continuous solution may instead be accompanied by a number of additional "weak discontinuity" solutions. Namely, solutions of the same boundary value problem that may represent fibre-scale failure modes.
It is recalled in this context that a preliminary theoretical analysis that enables determination of weak discontinuity surfaces in linearly elastic structural components reinforced by fibres resistant in bending is already available in [8,25]. Reference [8], in particular, makes use of precisely the same, restricted theoretical background employed in the present study, but applies its findings to three-dimensional study of composites reinforced by two families of straight uniaxial fibres. Specialisation of that analysis [8] to the case of a single fibre family is currently under way. This is expected to lead to weak discontinuity solutions that accompany potential three-dimensional extensions of the continuous solution detailed in Section 3. Identification of weak discontinuity solutions associated to the present, plane strain solution (Section 3) may then follow as particular cases.
where 1 2 ,  ,…, 5  are considered as known dimensionless constants, and a subscript "L" or "T" indicates the axis or the plane of transverse isotropy, respectively. With use of the mixture law (1), the effective elastic moduli of the functionally graded fibrous composite are obtained in the following form: 1 , while it is still (e.g., [26]) The elastic moduli appearing in (4) can then be obtained by inserting (A.3) and (A.4) into the standard relevant formulas (e.g., [26]), after aligning the longitudinal direction of transverse isotropy, L, with the x-axis of the adopted Cartesian co-ordinate system.

Appendix B: Implementation of the fictitious layers method
For a sufficiently thin plate or beam (h/L << 1), an approximate solution of (23) is obtained by replacing the variable z appearing in (6) and, hence, in (22) with its middle-plane value, namely its value at z = 0. In this manner, (23) is replaced by the following system of approximate linear ODEs: As this ODE system has constant coefficients, the exact form of its general solution can be expressed as follows: where the elements of the appearing exponential matrix S(z) are determined in the manner detailed in [19]. It is thus anticipated that the thinner is the inhomogeneous structural component of interest the nearer (B.2) approximates the exact solution of (23) or, equivalently, (20). When the thickness is not sufficiently small, the exact solution of (20) is approached computationally very closely by dividing the structure into N successive fictitious layers (see Figure 1) having the same constant thickness, h (j) = h/N (j = 1, 2,…, N). Each individual fictitious layer is associated with a local transverse co-ordinate parameter, z (j) = z -(j-1)h/N + h/2, and, due to the FGM nature of the plate, is itself materially inhomogeneous in the region - However, by choosing a suitably large value of N, each fictitious layer is itself regarded as a sufficiently thin plate or beam whose mechanical response and behaviour are described satisfactorily by an approximate solution of the form (B.2). The approximate solutions thus obtained for all N fictitious layers are then suitably connected together by means of appropriate continuity conditions imposed on the displacement and interlayer stress components. Upon increasing the value of N, this process provides a sufficiently close solution to that of the exact governing equations (20) (see also [18]).
In more detail, the continuity conditions imposed on a generic j-th material interface (denoted by zj in Figure 1) are as follows (j = 1, 2,…, N -1): In matrix form, these are transformed into the following:  N). Application of the same notation is extended to include the appearing fibre bending stiffness parameter, ( ) , where, however, it is also implied that previous use is made of (10). Upon using successively Eqs. (B.1), (B.2) and (B.4), one builds up the solution of the problem considered in a recursive manner. Hence, for the i-th layer, it is The value of ( ) on the outer lateral surface is then obtained as follows: If this is connected with the lateral boundary conditions (15), then (B.6) yields a linear system of four simultaneous algebraic equations for the four unknown components of the vector (−ℎ (1) /2) ≡ (−ℎ/2). Solution of that system of algebraic equations is then substituted back into (B.6) and provides a semi-analytical solution of the governing differential equations (23).
In the case of homogeneous fibrous composites, the first iteration of the outlined solution (N = 1) provides naturally the exact elasticity results obtained in [14]. For inhomogeneous composites, the number of iterations (N > 1) required for accurate prediction of displacement and stress distributions depends on the degree of the assumed material inhomogeneity. As is already mentioned, the convergence behaviour and success of this fictitious layer method has been verified repeatedly in the past (e.g,, [16][17][18][19][20]) as well as most recently in [21]. It accordingly suffices here to note that the value of  does not seem to exert significant influence on the observed convergence characteristics of the method, which thus remain essentially unchanged, regardless of whether the plate is made of polar ( 0   ) or non-polar material (λ = 0).
All numerical results shown in this communication were obtained by setting N = 100. In general, the maximum difference observed between corresponding results obtained on the basis of N = 100 and N = 70 iterations never exceeded 0.3%. It is worth noting that each iteration requires multiplication of 4×4 matrices only. As a result, the implied hundreds of iterations involved in computations do not require noticeable use of excessive computer time.

Appendix C: Consideration of the fibre-scale structure
It is assumed that fibres have circular cross-section of diameter d and, in the case of a homogeneous plate [14] are distributed along the z-direction in a regular form of N equidistant rows. The possible types of rectangular or triangular types RVEs depicted in Figure 3 consider that each vertex of an element is the centre of a fibre cross-section. In either case, Sy represents the distance of two neighbouring fibres in the y-direction.
Similarly, Sz represents the aforementioned constant distance between two neighbouring fibre rows. In this manner, Sz is the distance of neighbouring fibres in the zdirection of the rectangular element while, for a triangular RVE, it represents the height of the depicted isosceles triangle. In the particular case that the depicted triangle is considered equilateral (S = Sy), it is Sz = 3Sy/2.
It becomes then readily understood that, necessarily, the following conditions always hold: for the rectangular element. For the triangular element, these are modified as follows: For the purposes of the present study, d may be considered identical with the intrinsic parameter introduced in (9). However, the adopted notation distinction of those two parameters is retained here, in order to signify that (i) the shape of the fibre cross-section may be considered non-circular in different applications, and (ii) the intrinsic length parameter can acquire some different meaning in the theory of polar elasticity for fibre-reinforced materials [7], such as the fibre spacing for example.
Under these considerations, the fibre volume fraction of the RVE is defined as follows: For a rectangular RVE, this definition leads directly to This value of necessarily coincides with the maximum possible value of that the homogeneous counterpart of the present problem [14] is associated with when the fibre-scale structure is simulated with rectangular RVEs.
Similarly, maximum fibre volume fraction in a triangular RVE is achieved when S = Sy = d. In that case, (C.4) yields = 23 ≅ 0.907, (C. 6) which coincides with the maximum -value that the homogeneous version of the problem is associated with if the fibre structure is simulated with equilateral triangular elements. The values of shown in either (C.5) or (C.6) thus also consist corresponding upper limits that the -value can attain in the present inhomogeneous version of the problem, where fibres are assumed redistributed in the manner described by any of (31), (32) or (35) within the same matrix material. However, such an upper limit of can here be associated only with the densest fibre part of the inhomogeneous beam. Namely, the part located at the neighbourhood of z/h = 1/2 or z/h = -1/2 in top-stiff (31) or a bottom-stiff (32) beam, respectively, and the neighbourhood of both lateral planes (z/h = ±1/2) in the case of a beam reinforced in the symmetric manner (35). This is achievable by considering that Sz is a suitable function of z that takes its lowest value (namely Sz = d or Sz = 3d/2 for rectangular or triangular elements, respectively, in those densest fibre parts of the composite. Hence, by associating with the top (z/h = ½) or the bottom plane (z/h = -½) of a top-or bottom-stiff beam, respectively, either of (31) or (32) yields the maximum value of the parameter ε provided in (34). The inhomogeneous fibre distributions proposed in (31) and (32) are accordingly connected naturally with the present analysis when the fibre-scale structure is accurately simulated with rectangular or triangular RVEs, as soon as the εmaxvalue shown in (34 ) replaces the noted theoretical upper bound ε = 1. It is recalled that, by virtue of (3a), both (31) and (32) will thus still return < >=1/2. Table 1. Through thickness in-plane displacement distributions ̅ (0, ) of a top-stiff beam with volume fraction Vf = 0.5 + ε(z/h).