On three-dimensional dynamics of fibre-reinforced functionally graded plates when fibres resist bending

This communication aims to initiate an investigation towards understanding the influence that fibre bending stiffness has on the three-dimensional dynamic behaviour of fibrous composites with embedded functionally graded stiff fibres. In this context, it (i) formulates the general dynamical problem of a rectangular plate with embedded a single family of straight fibres that possess bending resistance and are distributed in a controlled, functionally graded manner through the plate thickness, and (ii) for simple support boundary conditions, it solves the free relevant vibration problem. The problem formulation is based on principles of polar linear elasticity and leads to a high-order set of Navier-type partial differential equations with variable coefficients. For simply supported edge boundaries, solution of these equations is achieved with the use of a computationally efficient semi-analytical (so-called fictitious layer) mathematical method. Two types of possible inhomogeneous distributions of straight fibres are considered for computational and numerical result presentation purposes. These are both regarded as possible, realistic types of inhomogeneous redistributions of stiff fibres that in previous studies have been assumed homogeneously distributed throughout the plate body. The presented numerical results examine to a considerable extent the manner that either of the employed types of inhomogeneous fibre redistribution, in conjunction with the fibre ability to resist bending, affects the dynamic behaviour of the fibrous composite plate of interest.

invention of carbon nanotubes and their use in composite structures as a type of fibre reinforcement gave rise to substantial interest on the influence that fibres exert on composite structures at a nanoscale order and dimension. This is because the low density and nanometer thickness of carbon nanotubes are coupled with remarkably high strength and stiffness and with similarly high bending resistance. However, initial mathematical modelling of fibrereinforced composites considered the fibres embedded in a structure as a kind of mathematical curves that possess neither thickness nor bending stiffness and resistance [1,2].
A theory that takes fibre thickness and relevant properties into consideration did not appear until relatively recently [3] . That theory dealt with modelling large and small elastic deformations of fibre-reinforced solids, and its linearized version has since redeveloped and extended independently [4,5] with the purpose to serve the study of the behaviour of structural fibrous composites with enhanced accuracy (e.g. [6][7][8][9]). By considering that there exist types of fibres that resist bending, the theory necessarily accounts for the effects of a couple-stress field and non-symmetric stress and thus is classified as a theory of polar linear elasticity. It is of interest to point out in this regard that Reference [6] identified a relatively simple stress analysis problem in which fibre existence can be accounted for only if the fibres embedded in the material of interest are considered capable to resist bending. In other words, conventional, non-polar linear elasticity fails in that case [6] to detect the presence of the fibres embedded in the material.
References [7,8] were thus devoted naturally to a number of relevant fundamental structural analysis problems, whose non-polar elasticity solution [10][11][12][13] had earlier provided substantial understanding and knowledge regarding both the static and the dynamic behaviour of fibrous structural composites. Those non-polar elasticity solutions [10][11][12][13] referred to the structural response of rectangular elastic plates containing fibres that are distributed either homogeneously or piece-wise homogeneously through the plate thickness. The corresponding non-polar elasticity solutions [7,8], achieved afterwards within the framework of the polar theory of linear elasticity of present interest [3][4][5], advance that knowledge and provide valuable information regarding the influence that very strong and stiff fibres, such as carbon nanotubes, exert in the behaviour of relevant fibrous composites.
Recognising the increasing interest in manufacturing and use of fibrous composites with functionally graded (FG) stiff fibres, Reference [9] initiated most recently an investigation towards understanding of the influence that fibre bending stiffness may have on the response of that type of inhomogeneous structural components. It is known that mathematical modelling of the controlled material inhomogeneity associated with FG composites leads to Navier-type elasticity equations that possess variable coefficients and, as a result, increases considerably the difficulty of their solution. Reference [9] thus started the implied investigation with an initial step that considered and investigated a fundamental plane strain problem relevant to the static behaviour of FG fibre-reinforced plates when fibres resist bending. The successful solution to that FG plate problem [9] led to initial collection of a substantial amount of useful relevant information and, in addition, provided mathematical and computational experience that enables consideration and solution of relevant three-dimensional boundary value problems.
In dealing with a relevant, three-dimensional, dynamic structural analysis problem, the present communication thus formulates and, for simple support boundary conditions, solves the free vibration problem of rectangular plates with embedded a single family of straight fibres that are distributed in a FG manner through the plate thickness and resist bending. As in [9], formulation of this problem is based on the relevant restricted version of polar linear elasticity [3,4] that involves only a single elasticity modulus of fibre bending resistance. The problem formulation is thus presented in Sect. 2 in a manner that resembles its Reference [9] counterpart, and leads to a high-order set of Navier-type partial differential equations with variable coefficients.
For the case of rectangular plates subjected to simple support boundary conditions, Sect. 3 employs next a suitable semi-analytical mathematical method (e.g. [6,9,[14][15][16][17][18]) that enables solution of that high-order set of differential equations in a computationally efficient manner. Sect. 4 introduces two types of possible inhomogeneous redistribution of straight fibres that, in earlier studies (e.g. [7,8,[10][11][12][13]15]), have previously assumed homogeneously distributed throughout the plate body. These types of possible realistic fibre redistribution are afterwards employed, in Sect. 5, for computational purposes that underpin the presentation and discussion of corresponding numerical results. That section thus examines to a considerable extent the manner that each of the employed types of inhomogeneous fibre distribution, in conjunction with the fact that fibres can resist bending, affects the dynamic behaviour of the

Theoretical formulation
A linearly elastic fibre-reinforced plate ( Fig. 1) has finite lengths, L 1 and L 2 , in the x and y directions, respectively, and finite thickness, h, in the z direction of a Cartesian coordinate system Oxyz 0 ≤ 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 through-thickness inhomogeneous distribution of the fibres is regulated by controlling their volume fraction, V f (z) which requires from the material properties of this structural component to be known functions of the z coordinate parameter. Every material property, P(z) say, of such a functionally graded fibrous composite can be defined according to 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. The inequality conditions noted in (1) are imposed on the ground of evident theoretical arguments that hold true regardless of the particular form of V f (z) or, equivalently,V m (z). In this context, the denoted upper limit of the fibre volume fraction, namely V f (z) = 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 they are distributed very densely.
The average fibre and matrix concentrations of the composite are defined as follows: The special case of a homogeneous fibrous composite, where the fibre volume fraction is constant, is thus characterised by the relationship V f (z) = V f for all z. If the fibres resist bending and V f is adequately high, say 40% to 60%, then fibre deformation generates considerable couple-stress and non-symmetric stress (e.g. [6][7][8][9]). In the case of FG fibrous composites with relatively low V f , creation of a couple-stress field is still possible locally, namely in specific parts of the composite where V f (z) anticipates high fibre concentration. Figure 2 illustrates the manner that the non-zero components of the deviatoric couple-stress tensor, along with the components of the non-symmetric stress tensor, act on an infinitesimal cuboid element of the elastic plate material.
The symmetric part of the stress tensor is given by the standard form of the generalised Hooke's law where the assumed transverse isotropy requires C 44 = (C 33 − C 23 )/2, and the linear elasticity strain components are Here, U (x, y, z, t), V (x, y, z, t) and W (x, y, z, t) are the displacement components along the x, y and z direction, respectively, and a comma denotes partial differentiation with respect to the indicated coordinate parameter(s). The elastic moduli appearing in (4) vary in the transverse direction in accordance with the mixture law (1), namely where, like in (1) as well as in what follows, the superfix m or f denotes parameters associated with the matrix or fibre phase of the composite, respectively. 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. The antisymmetric part of the stress tensor contains only two non-zero stress components. These are defined as follows [3][4][5]: where the appearing couple-stress components are and K f z and K f y represent the non-zero components of the fibre-curvature vector. Unlike C i j which have dimensions of stress, the fibre bending modulus d f has dimensions of force.
Like C i j , this fibre bending stiffness could also be expected to obey the mixture law (1). However, unlike the fibre phase, the isotropic matrix phase does not influence the fibre bending stiffness of the composite, and, as a result, the second term in the right-hand side of the corresponding expression (1.a) is zero. Hence, this modulus attains here the through thickness variable form where, in line with previous relevant studies [ [6][7][8][9], the intrinsic material length parameter l relates to the fibre thickness. It follows that l = 0 corresponds to the case of non-polar material behaviour, where fibres are perfectly flexible and the absence of couple-stress (m xy,x = τ [xz] = 0, m xz,x = τ [xy] = 0) enables the stress tensor to attain its conventional symmetric form.
In the absence of body forces, the equations of motion of three-dimensional polar elasticity acquire the form where ρ is the material density of the plate and a dot denotes differentiation with respect to time. Nevertheless, these equations are slightly simplified in the present case, where the antisymmetric part of τ yz is zero (τ [yz] = 0). With appropriate use of equations (4)- (11), these equations of motion attain a relevant Navier-type form, which is as follows: For analytical convenience, these partial differential equations are next rearranged in the following matrix form: where A = ⎡ ⎣ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ⎤ ⎦ , and the components of the matrix A are given as follows: The terms that appear underlined in expressions (15) are due to advanced modelling features that are not accounted for in conventional, non-polar elasticity studies that consider through thickness homogeneous fibre distribution (e.g. [10][11][12][13]15,17]). Terms underlined by a single line occur in Navier-type equations of non-polar linear elasticity only if fibres are still considered perfectly flexible but, also, distributed in a through-thickness FG manner. In a similar manner, terms underlined by a double or a triple line are all stemming from modelling features associated with effects of fibre thickness and fibre bending resistance, though the term underlined by the triple line attains a zero value if fibres are distributed homogeneously (e.g. [7,8]).

Free vibrations of a simply supported plate
One of the simplest applications that the outlined theoretical formulation can be related to is the free vibration problem of simply supported plates, namely elastic plates subjected to the following set of homogeneous boundary conditions at the edges x = 0, and x = L 1 : and the following set of homogeneous boundary conditions at y = 0 and y = L 2 : Due to the free vibration nature of the problem of interest, the lateral boundaries of the plate are assumed free of external tractions, thus leading to the following final set of homogeneous boundary conditions at z = ±h/2: It is pointed out that the number of the outlined boundary conditions is consistent with the order of the set of the Navier-type partial differential equations (12), which is ten in the x variable and six in either of the variables y and z.
A relevant comment is now added about the spherical part of the couple-stress tensor, which is customarily considered undetermined in this kind of Cosserat-type elasticity boundary value problems. It is accordingly assumed that, if the distribution of the spherical part of the couple-stress tensor could become available, through the solution of some relevant boundary value problem (e.g. [19]), then that distribution should enable each of the normal components of the couple-stress tensor, namely m 11 , m 22 and m 33 , to be zero value on any of the plate boundaries it acts upon.
Under these considerations, it is also noted that in the special case that fibre bending stiffness is neglected, the present problem reduces naturally to its non-polar elasticity counterpart (e.g. [20]). In further special cases that V f and V m are both known constants and, hence, the material of the fibrous composite plate is homogeneous, this problem reduces to its polar elasticity counterpart studied in [7]. As is already mentioned, corresponding non-polar elasticity problems that refer to homogeneous simply supported plates have already been considered and studied in [12,13].

Simplified form of the governing differential equations
It can readily be verified that the simple support boundary conditions (16) and (17) are satisfied exactly by the following choice of a displacement field: where the appearing natural frequency of vibration, ω, and the functions f (z) , r (z) and g (z) are to be determined. Expressions (19) thus represent plate vibration modes, each one of which is a potential solution to the described boundary value problem. Upon inserting (19) into (13), the latter equation is transformed into a set of simultaneous ordinary differential equations (ODEs) with variable coefficients. This can be expressed as follows: where The remaining of the entries appearing in (21a) are generally functions of z and are given in Appendix A.

Solution with use of the fictitious layer method
Solution of (20) is next achieved with the use of a semi-analytical method, which considers that the inhomogeneous polar material plate 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 Fig. 1), the larger the number of those fictitious layers, the nearer the obtained numerical results approach their exact elasticity counterparts. This method has been applied successfully to relevant static and dynamic studies of homogeneous and laminated composite components (e.g. [6,9,[14][15][16][17][18] and relevant references therein), and its numerical stability and the rate of convergence are found superior to those of corresponding analytical solutions based on power-series methods.
Description of the solution thus obtained is facilitated by initially converting (21) into the following equivalent set of six first-order linear ODEs with variable coefficients: where and the elements of the matrix T are given in Appendix A. For a sufficiently thin plate (h/L 1 << 1), an approximate solution of (22) is obtained by replacing the variable z appearing in (6) and, hence, in (21) with its middle-plane value, namely z = 0. In this manner, (22) 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 [15][16][17][18]. The thinner is the inhomogeneous plate of interest the better (25) thus approximates the exact solution of (22). When the thickness is not sufficiently small, the exact solution of (22) is approached computationally very closely by dividing the structure into N L successive fictitious layers (Fig. 1) having the same constant thickness, Each individual fictitious layer is associated with a local transverse coordinate parameter, z ( j) = z − ( j−1)h N L + h 2 , and, due to the FGM nature of the plate, is itself materially inhomogeneous in However, by choosing a suitably large value of N L , 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 (25). The approximate solutions thus obtained for all N L fictitious layers are then suitably connected together by means of appropriate continuity conditions imposed on the displacement and the interlayer stress components. Upon increasing the value of N L , this process provides a sufficiently close solution to that of the exact governing equations (20) (see also [16]).
The continuity and interface equilibrium conditions imposed on a generic j-th material interface (denoted by z j in Fig. 1) are then as follows ( j = 1, 2, . . . , N L − 1): In matrix form, these are transformed into the following: where and 13 , …, etc, signify the constant values that the implied elastic moduli acquire on the middle plane, z ( j) =0, of the j-th fictitious layer ( j = 1, 2, . . . , N L ). Application of the same notation is here extended to include the appearing fibre bending stiffness parameter, d f ( j) , where, however, it is also implied that previous use is made of (10).
Upon using successively Eqs. (26)-(28), one builds up the solution of the problem considered in a recursive manner. Hence, for the i-th layer, it is where The value of F (z) on the outer lateral surface is then obtained as follows: By considering the boundary conditions (18), eqs. (31) yield a system of 6N L homogeneous algebraic equations for the 6N L unknowns, C . This can be written in the following form of a generalised eigenvalue problem: where, as is implied by (A.4) and (25), the components of the matrix Γ are given in Appendix A and are generally dependent on the square of the frequency parameter, ω. Solution of this generalised eigenvalue problem requires from the determinant of the matrix Γ to be set equal to zero. The resulting transcendental characteristic equation has an infinite number of roots. Any real and positive of those roots represents the square of a natural vibration frequency. Following its numerical determination, a subsequent substitution of the value of any of those frequencies back into (29) enables determination of the corresponding eigenvector, which represents the relevant mode shape of the freely vibrating plate.

Application for selected forms of the fibre volume fraction
As is mentioned in Sect. 2, and noted also in (2), the maximum fibre volume fraction, V f max , is always less than 1 in real engineering applications. As pragmatic values of V f max necessarily depend on the fibre structural architecture, Reference [9] provides a detailed discussion of the manner that this can be evaluated with the use of appropriate representative volume elements (RVE) of the fibre distribution pattern [21]. With the help of Fig. 3, that discussion, which is not repeated here, revealed that On the other hand, potential applications of the present analysis are facilitated through its plausible connection with relevant studies of corresponding homogeneous fibre-reinforced plates. It is accordingly recalled [9] that the material properties employed in relevant numerical applications, performed in [7] for corresponding homogeneous fibre-reinforced plates, are as follows: It is thus fitting to mention at this point that, for simply supported plates with homogeneously distributed fibres, the present analysis gives identical natural frequencies with those obtained in [7]. Moreover, as in that case the plate is homogenous rather than a functionally graded one, those results are obtained with the first iteration of the outlined computational process (N L = 1). It is noted in this regard that the outlined verification of the present analysis is possible for any constant choice of V f =< V f >, as long as the individual elastic moduli of the fibre and matrix phases enable (34) to represent the effective material properties of the fibrous composite plate. In accordance with the principal interest of the present analysis, this Section thus introduces next two types of a possible inhomogeneous redistribution of the straight fibres that (i) were assumed homogeneously distributed in earlier studies (e.g. [7,8,[10][11][12][13]15]), and (ii) their homogenous through thickness distribution enabled (34) to represent the effective elastic moduli of the fibrous composite plate.

Linear redistribution of the fibres
To study the effect of nonuniform fibre distributions on the plate vibration frequencies, the following linearly inhomogeneous fibre distributions are initially considered (see also [9]): In this case, (3a) yields and is thus seen that (35) represents an inhomogeneous composite consisting of 50% fibre and 50% matrix material. The fibres are distributed linearly along the thickness direction and top or bottom stiff plates are obtained according to ε > 0 or ε < 0, respectively. It should be noted that the special case ε = 0 represents homogeneous fibre distribution considered in [7]. The corresponding fibre distributions are depicted in Fig. 4. However, the stress analysis detailed [9] is affected by the manner that the external loading is applied on the plate boundaries, while in the present, free vibration problem the plate is not subjected to any kind of externally applied loading. It follows that the free vibration problem of a top stiff and a bottom stiff plate resulting from the In principle, ε max may be as high as 1 but, in real engineering applications, the fibre-scale structure of a composite relates this parameter with V f max . Connection of (35) with (33) thus reveals that maximum volume fraction is obtained at the top V f h 2 = 1+ε 2 and minimum at the bottom V f − h 2 = (1 − ε)/2 lateral surface of the plate. Hence, as is also detailed in [9], it is found that .5708 for rectangular RVEs, 0.812 for triangular RVEs.
The relevant numerical results presented in Sect. 5 below refer to inhomogeneous composites whose effective material properties are those detailed in (34). These are obtained by connecting the mixture law (1) with the linear fibre distribution (35) and considering that the following relationships hold between the elastic moduli of the where the superscript " f " indicates elastic moduli associated with a transversely isotropic fibre phase, while E and ν stand for the Young modulus and the Poisson ratio, respectively, of the corresponding isotropic matrix phase. More details of the inverse process in which (38) leads to (35) are provided in [9]. 4.2 Through-thickness symmetric, piece-wise linear redistribution of the fibres Symmetric distributions of fibres with respect to midplane of the plate are also possible [9], and a piece-wise linear such can be obtained through use of the following expression of a fibre volume fraction: where α is some real positive constant. In this type of fibre distribution, the fibre volume fraction is zero at the midplane and, while it increases linearly in moving towards either of the lateral plane directions, it reaches its maximum value at those top and bottom planes V f h 2 = V f − h 2 = α 2 . However, as is also shown in [9], α = 2 is not possible for either rectangular or triangular fibre RVE arrangements. Instead, the maximum volume fractions noted in (33) are reached for the following maximum value of α: verify that neither of the values of α noted in (40) enables consideration of an average fibre volume fraction which is as high as its 0.5 counterpart assumed with the linear fibre distributions (35). However, as is detailed in [9], the effective material properties (34) can still be obtained for an average volume fraction and V f = 0.45, by replacing (38) with the following: The numerical results presented in the next section for symmetric fibre distributions are thus obtained by using the effective material properties (34), and α = 1.814 in (39). This is the maximum value that α can obtain in the case of the triangular RVE architecture.

Numerical results and discussion
As is already mentioned, by setting ε = 0 in (35), one obtains the homogeneous fibre distribution considered in [7]. All relevant numerical results obtained with the use of the present analysis for ε = 0 were found in excellent agreement with their counterparts tabulated in [7]. In this context, material parameters that refer to composite plates having their embedded fibres uniformly distributed are distinguished by the suffix "uni". All numerical results obtained for ε = 0 thus refer to elastic plates whose functionally graded inhomogeneity is formed through a specified redistribution of the same volume of fibres encountered in a corresponding plate with uniformly distributed fibres.
The numerical results presentation is facilitated through use of the following non-dimensional frequency parameter: As is already mentioned, for each pair (m, n) of longitudinal and transverse half-wave numbers, the roots of the characteristic equation stemming from (29) yield an infinite number of frequencies. The results presented in what follows are restricted to specific, non-negative integer values of m and n that do not exceed 2.
The discussion thus focusses on natural vibration modes whose natural frequency parameter conveniently represented by the notation Ω 01 , Ω 10 , Ω 11 , Ω 21 , Ω 12 and Ω 22 . Moreover, from the essentially infinite number of frequencies that may be determined for each of the implied pairs of (m, n)-values, attention focuses only on the lowest one, which is the most likely frequency to cause resonance failure of the composite structural component.
It is also noted that, as expressions (19) make it clear, in either of its Ω 01 and Ω 10 modes, the plate vibrates in a purely in-plane, one-dimensional manner. This is because, in either of these modes, the transverse displacement component, W , and one of the in-plane displacement components, namely V or U , respectively, is zero. Obviously, there exists no Ω 00 vibration mode as, in that case, all three displacement components are equal to zero.
In line with previous relevant studies [8,9], the following non-dimensional parameter: is associated with the intrinsic material parameter, l, that refers to fibre thickness. As is also noted in [9], this parameter should not be misinterpreted as denoting the Lamé modulus employed in (7). Further discussion regarding the nature of the non-dimensional fibre thickness parameter (44), as well as of possible restrictions that can be imposed of its values, is given in [9] and will not be repeated here. However, it is appropriate to mention that by virtue of (10), (44) essentially results in the following reparametrisation the fibre bending stiffness modulus: and that, naturally, λ = 0 corresponds to the special case of non-polar material behaviour where fibres are considered perfectly flexible.

Convergence of the fictitious layer method
As an example of the convergence features of the fictitious layer method employed in Sect. 3.2, Table 1 demonstrates the manner that the value of the aforementioned dimensionless frequency parameters (43) of a square, top-stiff functionally graded fibre-reinforced plate (ε = 0.1) varies with increasing the value, N L , of fictitious layers. A comparison between corresponding numerical results obtained with the use of seven and twenty-one fictitious layers shows that, regardless of the value of the non-dimensional fibre thickness parameter, λ, consideration of seven layers suffices for the chosen plate thickness (h/L 1 = 0.1) to provide numerical results which are accurate in at least three significant figures. For much thicker plates, or for plates with considerably higher inhomogeneity, achievement of the same level of accuracy may require an increased value of the N L number. Nevertheless, all numerical results shown next are obtained for plates with the same plate thickness to fibre-length ratio (h/L 1 = 0.1) and, in all cases, the use of twenty-one fictitious layers (N L = 21) was found adequate to provide results which are very accurate for practical purposes.
It was also observed that the convergence of the fictitious layer method was substantially faster in the present dynamic analysis case than it was in the case of the stress analysis considered in [9]. This difference is attributed to the well-known fact that accurate prediction of global response characteristics of thin-walled structures, such as frequencies, buckling loads and through-thickness averaged displacements, is always easier achievable than it is for corresponding stress, strain and displacement distributions. The latter are strongly influenced by the local structural response and are thus also classified as local or detailed response structural characteristics of the thin-walled structure of interest.

Vibration frequencies for top-stiff fibre-reinforced plates
In this context, Table 2 presents the convergent values of the lowest frequency parameters, Ω mn , obtained for the same square plate considered in Table 1, along with the corresponding frequency values obtained when that plate is either homogeneous (ε = 0) or its inhomogeneity attains the highest value possible for triangular RVEs, namely ε = 0.812; see (37). Corresponding numerical results for rectangular plates with in-plane aspect ratio L 2 /L 1 = 1/2 and L 2 /L 1 = 2 are tabulated in Tables 3 and 4, respectively. It is recalled in this context that, in all cases, the fibres are aligned with the longitudinal, x 1 -direction, and L 1 thus represents their length as well.
The remaining of the frequency parameters, where m > 0 and n > 0, are all associated with vibrational modes that couple completely both in-plane and out-of-plane motions; see (19). It is observed that all these frequencies are increasing with increasing the fibre bending stiffness parameter, λ. The relative difference between corresponding frequencies predicted with the use of conventional (non-polar) elasticity and its present, polar elasticity extension attains its highest increase (approximately 100% for λ = 0.01) for m = 2 and is getting higher by increasing L 2 /L 1 ratio. This is the case in which the fibre bending motion involves two half-waves and, hence, meets naturally the highest fibre bending resistance. On the other hand, the observed effect of fibre bending stiffness appears to be decreasing with decreasing the L 2 /L 1 ratio. The lowest fibre bending stiffness effect is observed in the case that L 2 /L 1 = 0.5.
In the case of square plates (Tables 1 and 2) the frequency Ω 21 is always higher than its Ω 12 counterpart regardless of whether the fibres are perfectly flexible or resistant in bending. This is because, despite the implied geometrical symmetry in the x-and y-directions, the presence (absence) of strong fibres in the longitudinal (transverse) direction increases (decreases) resistance of half-wave creation in the transverse (longitudinal) direction. This observation is not affected by changes of the inhomogeneity parameter, ε, and still holds for the frequencies of rectangular plates with L 2 > L 1 (e.g. Table 4). However, as is seen in Table 3, where L 2 = L 1 /2, this is not always the case if L 2 < L 1 . The results presented in that Table provide a clear indication that, for L 2 < L 1 , the order of the frequencies Ω 21 and Ω 12 may well be affected by the inhomogeneity of the fibre distribution, or the fibre bending stiffness, or a combination of these fibre-related material properties. For instance, Ω 12 in Table 3 is higher than Ω 21 in the non-polar elasticity case (λ = 0), while the opposite is observed when λ is non-zero.
No further general trends appear easily observable for the top-stiff plate material configuration through the results presented in Tables 1, 2, 3, and 4. Nevertheless, Table 5 presents the corresponding vibration frequency parameters obtained for a top-stiff functionally graded rectangular plate in a limiting case that the plate geometrical characteristics enable resemblance of corresponding plane strain deformation. It is recalled that Reference [9] is already devoted to a corresponding, static, stress analysis problem.
It is observed that, with the exception of Ω 01 , the frequency parameters shown in Table 5 do not substantially deviate from their counterparts tabulated in Table 2 for L 2 /L 1 = 2. It is thus concluded that, due to the strong fibre reinforcement that the plate is equipped with in its longitudinal direction, it approaches its transverse plane strain limit not far beyond L 2 = 2L 1 As an exception to this observation, Ω 01 appears to be continuously decreasing with increasing L 2 /L 1 but, as is explained in Appendix B, this is only due to the effect that the parametrisation (43) has on the particularly simple solution that the problem meets in that case

Through-thickness symmetric fibre redistribution
In dealing with the through-thickness symmetric, piece-wise linear fibre redistribution considered in Sect. 4.2, the present analysis employed further the form (39) of fibre volume fraction in conjunction with the maximum value that the real positive constant α attains when fibres are arranged in accordance with the triangular RVE specification; see (40). In this context Table 6 presents relevant numerical results that correspond to their top-stiff plate counterparts tabulated in Tables 2, 3, 4, and 5. In more detail, the second part of Table 6 presents non-dimensional frequencies that correspond to the square plate frequency parameters tabulated in Table 2. In a similar manner, the results shown in the top and the third the parts of that Table correspond to their counterparts presented in Tables 3 and 4, respectively, while the contents of the last part of Table 6 match their practically plane strain counterparts demonstrated in Table 5.
The numerical results presented in Table 6 show that, for the reasons explained in Appendix B, the frequency parameters Ω 01 and Ω 10 behave in a manner similar to that observed for their counterparts tabulated in Tables 2, 3, 4, and 5. However, the remaining frequency parameters show different trends to those observed previously for their top stiff plate counterparts. In fact, unless the plate is wide enough (L 2 /L 1 > 2), even the frequency parameter Ω 21 seems to be smaller than its Ω 12 counterpart in most cases. In this regard, it is quite remarkable that, in the limiting case of plane strain, most values of the frequency parameters tabulated in the bottom part of Table 6 are comparable with their top-stiff plate counterparts shown in Table 5.

Conclusions
The numerical results presented in the preceding section naturally confirm the expectation that a different redistribution of initially homogeneously distributed fibres within the same matrix material of an elastic plate can considerably affect the dynamic characteristics of the plate. This fact holds of course regardless of whether the fibres embedded in a fibrous composite plate are perfectly flexible or resistant in bending. Moreover, it naturally holds true not only in cases that, like the present one, refer to plate dynamics, but also to relevant static, stress analysis cases, such as the one considered in Reference [9]. It is emphasised in this context, that the fictitious layer method employed here as well as in [9] enables fast converging solution of relevant boundary value problems, regardless of the form of the variable coefficients appearing in the governing differential equations. In this context, further types of fibre distribution, in addition to the few ones employed already within the same or some different rectangular simply supported composite plate, can still be studied with relatively ease, through the use of the three-dimensional polar elasticity analysis detailed both in [9] and in this communication. It should not be overlooked though that plates with all around simply supported edges are mathematically privileged structures, in the sense that they enable the Navier-type governing equations to free themselves form in-plane coordinate dependence in a straightforward manner. For more realistic types of boundary conditions, three-dimensional dynamic or static analyses of the functionally graded fibrous composites of present interest become considerably more demanding. Consideration of different sets of edge boundary conditions becomes, in fact, particularly more challenging within the here considered framework of anisotropic polar linear elasticity that anticipates emergence and influence of couple-stress and non-symmetric stress fields. It is recalled in this connection that, within the conventional framework of non-polar linear elasticity, threedimensional boundary value problems of plates subjected to different types of edge boundary conditions are alternatively tackled with the use of the energy minimisation methods, such as the method of Ritz. While relevant energy minimisation theorems, as well as the manner they are used in non-polar elasticity applications, are available for long time, such theorems started to appear only recently in relation to polar linear elasticity of fibrous composites (see [19] and references therein). These new developments are expected to assist formation of appropriate mathematical tools that will enable extension of the analysis presented here as well as in [9] towards the study and understanding of the mechanics of polar fibre-reinforced plates subjected to different sets of edge boundary conditions. Still though, the simply supported plate solutions presented here as well as in [7] and [9] are particularly useful when they are suitably combined with appropriately developed two-dimensional, polar elasticity models of thin-walled structural components (e.g. [22]). It is recalled in this connection that, by incorporating a non-polar elasticity solution for simply supported plates into the transverse deformation shape functions of a corresponding two-dimensional elasticity model, the latter is enabled to return (i) exact elasticity predictions for simply supported plates [23,24], and (ii) very accurate elasticity predictions for either simply supported thin shells or thin plates and shells subjected to different sets of edge boundary conditions [25,26].
In dealing with polar linear elasticity of fibrous composite plates, such a combination of a two-dimensional plate model [22] with a corresponding three-dimensional simply supported solution is already applied successfully in cases that fibres resistant in bending are homogeneously distributed through the plate thickness [7]. The simply supported plate solutions presented in [9] as well as in the present communication can then play the same valuable role in the case that thin polar structural components reinforced with functionally graded fibres resistant in bending. Moreover, by neglecting the effect of fibre bending stiffness, that combination will naturally attain a simpler form that is still similarly valuable in the mechanics of functionally graded thin-walled structures reinforced with perfectly flexible fibres. In dealing with the in-plane distortional mode (m, n) = (0, 1), the displacement field (19) yields and this enables the equations of motion (12) or, equivalently (14) to reduce to the single equation Due to the plate inhomogeneity, (B.2) is generally a second-order ordinary differential equation with variable coefficients that can be solved either with standard power-series methods or with the fictitious layer method employed in this communication.
However, in the case of a homogeneous plate, where C 55 is constant, (B.2) simplifies further and becomes The lowest frequency associated with this mode is thus seen associated with the linear mode shape where A is a constant. This solution of (B.2) returns the natural frequency ω = π L 2 C 55 ρ , (B.4) which, by virtue of (43), obtains the non-dimensional form The fact that all Ω 01 -values illustrated in Tables 1, 2, 3, 4, and 5 for homogeneous plates (ε = 0) can alternatively be obtained with direct use of (B.5) verifies the efficiency and correctness of the employed computational code. The latter was naturally used for the evaluation of the remaining of the Ω 01 -values shown in Tables 1, 2, 3, 4, and 5 for ε = 0. These values clearly demonstrate a small influence that the assumed, top-stiff plate inhomogeneity exerts on the values (B.5) of their homogeneous plate counterparts. Nevertheless, Table 6 reveals that the Ω 01 -values are influenced more severely from the enhanced inhomogeneity encountered in a through-thickness symmetric fibre redistribution.
Similar observations apply with regard to the in-plane distortional mode (m, n) = (1, 0), for which the displacement field (19)  which, due to the plate inhomogeneity, is generally again a second-order ordinary differential equation with variable coefficients. It should be noted though that, unlike (B.2), (B.7) is now influenced by resistance that the fibres may exhibit if/when subjected to in-plane bending. However, in the case of a homogeneous plate, where C 44 is constant, (B.7) simplifies further and becomes The lowest frequency associated with this mode is again associated with the linear mode shape where B is a constant. This solution of (B.8) returns the natural frequency and, by virtue of (43), obtains the non-dimensional form (B.11)