Static inconsistencies in certain axiomatic higher-order shear deformation theories for beams, plates and shells.

Static inconsistencies that arise when modelling the ﬂexural behaviour of beams, plates and shells with clamped boundary conditions using a certain class of axiomatic, higher-order shear deformation theory are discussed. The inconsistencies pertain to displacement-based theories that enforce conditions of van- ishing shear strain at the top and bottom surfaces a priori . First it is shown that the essential boundary condition of vanishing Kirchhoff rotation perpendicular to an edge ( w ; x ¼ 0 or w ; y ¼ 0) is physically inac- curate, as the rotation at a clamped edge may in fact be non-zero due to the presence of transverse shear rotation. As a result, the shear force derived from constitutive equations erroneously vanishes at a clamped edge. In effect, this boundary condition overconstrains the structure leading to underpredictions in transverse bending deﬂection and overpredictions of axial stresses compared to high-ﬁdelity 3D ﬁnite element solutions for thick and highly orthotropic plates. Generalised higher-order theories written in the form of a power series, as in Carrera’s Uniﬁed Formulation, do not produce this inconsistency. It is shown that the condition of vanishing shear tractions at the top and bottom surfaces need not be applied a priori , as the transverse shear strains inherently vanish if the order of the theory is sufﬁcient to capture all higher-order effects. Finally, the transverse deﬂection of the generalised higher-order theories is expanded in a power series of a non-dimensional parameter and used to derive a material and geometry dependent shear correction factor that provides more accurate solutions of bending deﬂection than the classical value of 5/6. (cid:2) 2014 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Literature Review
In practical applications composite laminates are typically modelled as thin plates and shells because the thickness dimension t is an order of magnitude smaller than representative in-plane dimensions l x and l y . This feature allows the problem to be reduced from a three-dimensional (3D) to a two-dimensional (2D) one coincident with a chosen reference surface of the plate or shell. The major advantage of this approximation is a significant reduction in the total number of variables and computational effort required. These theories are aptly called Equivalent Single Layer (ESL) theories as the through-thickness properties are compressed onto an equivalent single layer by integrating properties of interest through the thickness. Many ESL theories are based on the axiomatic approach, whereby an intuitive postulation of the displacement and/or stress fields in the thickness z-direction is made. Appropriate displacement-based, stress-based or mixed variational formulations are then used to derive variationally consistent governing field equations and boundary conditions.
The most prominent example of the ESL theories is the Classical Theory of Plates (CTP) developed by Kirchhoff [1] and revisited by Love [2]. The principle assumptions of Kirchhoff's theory are that: where the comma notation signifies differentiation. This approach reduces the number of displacement fields to three, namely the two membrane displacements u 0x and u 0y and the constant transverse deflection w of the neutral plane.
Kirchhoff's theory is currently the most widely used formulation for analysing thin-plate structures and also forms the foundation for Classical Laminate Analysis (CLA) [3]. In reality, the analysis of layered composite structures is significantly more cumbersome due to a plethora of new features such as in-plane anisotropy (IA); transverse anisotropy (TA); and interlaminar displacement, transverse shear and transverse normal stress continuity (IC). For example, composite laminae often exhibit much greater values of Young's modulus orthotropic ratio (E 11 =E 22 ¼ E 11 =E 33 % 140=10 ¼ 14), i.e. in-plane orthotropy, than isotropic materials (E xx ¼ E yy ¼ E zz ). Furthermore, the induced error of Kirchhoff's hypothesis for isotropic plates is approximately 5% for thickness to characteristic length ratios of about 1/10 [4]. However, composite laminates exhibit much greater transverse orthotropy such that the longitudinal to transverse shear modulii ratios are approximately one order of magnitude greater than for isotropic materials (E iso =G iso ¼ 2:6; E 11 =G 13 % 140=5 ¼ 28). Thus, the transverse shear and normal stress deformations have a much greater effect on the flexural behaviour of the plate [4,5]. The nondimensional ratio k ¼ E 11 =G 13 ðt=LÞ 2 drives what Everstine & Pipkin [6] coined the ''stress-channelling'' effect on axial stress as shown in Fig. 1, where t is the thickness and L a characteristic planar dimension. As the value of k increases the cross-section shears increasingly and therefore transitions from a constant rotation to a higher-order distortion field. As a result, axial stress is channelled towards the outside surfaces.
Furthermore, transverse orthotropy, i.e. the difference in layerwise transverse shear and transverse normal modulii, leads to sudden changes in slope of the three displacement fields u x ; u y and u z at layer interfaces known as the Zig-Zag (ZZ) phenomenon. While interlaminar continuity of the displacements at layer interfaces requires u x ; u y and u z to be C 0 continuous at interfaces, interlaminar continuity of the transverse stresses leads to a C 1 discontinuous displacement field. In fact Demasi [8] showed that the ZZ form of the in-plane displacements u x and u y can be mathematically demonstrated from s xz and s yz , respectively while the ZZ form of u z can be derived from r z . Therefore an accurate model for multi-layered composite and sandwich structures should ideally address the modelling issues named C 0 z -requirements by Carrera [5,9]: 1. Through-thickness z-continuous displacements and transverse stresses (called interlaminar continuity, IC), 2. Discontinuous first derivatives of displacements between layers with different mechanical properties, i.e. the Zig-Zag effect, without becoming overly computationally expensive. In fact Carrera points out that ''compatibility and equilibrium, i.e., ZZ and IC, are strongly connected to each other'' [10]. However, the broad scope of work published in recent decades has shown that in many practical engineering cases, where TA is consciously minimised by lamination guidelines, the axiomatic theories developed for isotropic or single-layer structures are also applicable for laminated structures [11].
Since the first half of the 20th century a number of models have been presented that partially or completely revoke at least one of Kirchhoff's original assumptions. One of the earliest examples is the First-order Shear Deformation Theory (FSDT) which assumes that normals to the reference surface do not necessarily remain normal after deformation. To capture this the rotations h x and h y of the cross-section with respect to the undeformed state are introduced as new degrees of freedom, In FSDT the effect of shear deformation on the cross-section is captured in an average sense. Timoshenko [12] famously applied this hypothesis to the classical model for isotropic Euler-Bernoulli beams, while a two-dimensional extension for isotropic and singlelayer plates was presented by Mindlin [13] and extended to multilayered plates by Yang et al. [14]. FSDT represents an improvement on global structural phenomena such as bending displacement and low-frequency buckling and vibrational modes. However, FSDT does not improve localised strain and stress predictions, especially for highly heterogeneous and thick composite and sandwich laminates [15], because of its limiting assumption regarding uniform transverse shear strain. Furthermore, FSDT produces piecewiseconstant transverse shear stresses that violate continuity at layer interfaces and do not disappear at the top and bottom surfaces. Finally, as the actual transverse shear stress profile is at least quadratic in order shear correction factors are needed to energetically adjust the constant through-thickness strain profile of FSDT. Determining the magnitude of these shear correction factors is not a straightforward task, especially in the case of highly heterogeneous laminates, and various methods addressing such concerns have been published in the literature [16][17][18].
To take into account the actual higher-order distribution of introduced into the in-plane expansion u x thereby reducing the number of variables to that of Timoshenko beam theory. The associated bending moment and shear force were substituted into the well-known beam equilibrium equations, resulting in a variationally inconsistent fourth-order differential equation featuring only the transverse displacement w. Reddy [21] extended Levinson's theory to two-dimensional problems featuring u x ; u y and u z and derived the governing field and boundary equations in a variationally consistent manner using the principle of virtual displacements.
The two steps followed by Levinson [20] and Reddy [21] of firstly allowing the cross-section to distort in a higher-order fashion and secondly enforcing the physical boundary conditions of vanishing transverse shear strains at the top and bottom surfaces, has led to a class of theories that may collectively be written as Here the unknown variables c i capture the magnitude of the crosssectional distortion, where f ðzÞ is a pertinent shape function that approximates the parabolic distribution of the transverse shear strains through the thickness while guaranteeing that transverse shear strains vanish at the surfaces. A large number of different shear shape functions f ðzÞ have been published ranging from polynomial [22][23][24] to trigonometric [25][26][27][28][29], hyperbolic [30,31] and exponential [32,33] some of which are shown in Table 1. Karama et al. [32] pointed out that an expansion based on the combinations of an exponential and linear function in z is superior to a trigonometric function as it has all even and odd powers in the expansion and the coefficients of higher-order terms do not decay as quickly.
Another type of HOT known as the Refined Plate Theory (RPT) was introduced by Shimpi [34]. In RPT the transverse displacement is split into separate bending and shear components and the axial displacement is a function of individual bending and shearing rotations. As a result, the bending components do not contribute towards transverse shear forces and shearing components do not contribute towards bending moments. The advantages of this approach are that the governing equations maintain an intuitive resemblance to CPT and allow the development of finite elements that are free from shear locking.
Transverse normal strains may be incorporated by extending the expansion of the out-of-plane displacement u z to yield a class of theories denoted as Advanced Higher-Order Theories (AHOT). Here the class of theory is generally denoted by fa; bg where a refers to the order of expansion of the in-plane displacements u x and u y , and b to the order of the transverse displacement u z . These theories thus take account of Koiter's Recommendation which postulates that refinements of CPT need to account for both transverse shear and normal strains [35]. Milestone works were presented by Hildebrand, Reissner & Thomas [36] and Lo, Christensen & Wu [37], where the latter is a f3; 2g AHOT given by Further examples of AHOT are given by [37][38][39][40]. Generally, AHOT only provide improvements that are worth their additional computational effort for sandwich panels with compliant, thick cores [8] or when one face-laminate is considerably stiffer than the other [41]. Also note that Eq. (4) is a generalised expansion of the displacement field ðu x ; u y ; u z Þ in terms of the transverse variable z. Thus the boundary condition of vanishing transverse shear strains at the top and bottom surfaces is not enforced a priori and the rotations w ;i do not arise as in Eq. (3a).
The notion of a generalised axiomatic expansion was developed further by Carrera in what is known as the Carrera Unified ormulation [9,42,43] and its extension the Generalized Unified Formulation by Demasi [44]. Carrera's Unified Formulation is a hierarchical formulation where Taylor series or Lagrange polynomials are used to approximate the displacement fields over the cross-section [43]. This allows the order of the theory to be expressed as an input to the analysis. In this manner, the governing field equations are formulated based on the generalised axiomatic expansion. Theories of different orders can thus easily be applied without the need for separately deriving new field equations. The Taylor series or Lagrange polynomial expansion is not strictly limited to a function of the z-coordinate such that bi-axial bending, torsion and warping can be modelled by expanding in the x-and y-coordinates. Thus, in the framework of Carrera's Unified Formulation the displacement field is expressed as the expansion of generic functions F s , where F s are individual functions of the coordinates x, y and z;ũ s is the displacement vector; M is the number of terms in the expansion and according to the Einstein notation repeated indexes denote summation. For example, This representation allows the stiffness matrix of the theory to be expressed in terms of a few fundamental nuclei which are through-thickness integrals of material stiffness terms multiplied by a combination of F s terms.
All of the previously discussed theories are based on displacement formulations where the displacements u x ; u y and u z are treated as the unknown variables. Consequently, all strains and stresses are derived from the displacement assumptions using kinematic and constitutive equations, respectively. Governing field equations are typically derived in a variationally consistent manner using the Principle of Virtual Displacements (PVD). A disadvantage of these theories is that accurate transverse stresses are not easily derived from the constitutive equations [10].
Another class of models is based on applying the Hellinger-Reissner mixed variational principle. Here the strain energy is written in complementary form in terms of in-plane and transverse stresses, and the transverse equilibrium equation is introduced as a constraint condition using a Lagrange multiplier. Using this framework Reissner [45,46] developed an enhanced theory from the piecewise-linear in-plane stresses of CPT and parabolic transverse shear equations derived from Cauchy's 3D equilibrium equations. Recently, this approach has been applied to symmetrically laminated straight-fibre composites [47], variable stiffness composites [48] and sandwich panels where ZZ effects are important [49]. The advantage of these mixed variational statements is that both axial and transverse shear stresses can generally be captured accurately directly without the need for further post-processing steps. Reissner [50], considering multi-layered structures, had the insight that it is sufficient to restrict the stress assumptions to the transverse stresses because only these have to be specified independently to guarantee the IC requirements. This variational statement is known as Reissner's Mixed Variational Theory (RMVT) which makes model assumptions on the three displacements u x ; u y and u z and independent assumptions on the transverse stresses s xz ; s yz and r z . Compatibility of the transverse strains is enforced by means of Lagrange multipliers. Murakami [51] was among the first to apply RMVT to composite plates and, by introducing a piecewise continuous ZZ function that alternatively takes the values of +1 or À1 at layer interfaces, satisfied both IC and ZZ requirements. Recently, the RMVT has been applied with more physically accurate ZZ functions in the Refined Zigzag theory [41] introduced by Tessler et al. [52]. Extended reviews of ZZ theories and the application of RMVT to layered structures are provided by Carrera [10,5,53].

Contents of the paper
No matter the choice of shape function f ðzÞ, all theories written in the form of Eq. (3a) share the common feature that the CPT rotations w ;i feature in the axiomatic expansions of the displacement fields u x and u y . The presence of these terms leads to essential boundary conditions dw ;x ¼ 0 and dw ;y ¼ 0 and associated natural boundary conditions on the higher-order moments, when the governing field equations are derived in a variationally consistent manner using the principle of virtual displacements. Restraining the rotations perpendicular to clamped edges, ie. w ;x ¼ 0 and/or w ;y ¼ 0 as done by many authors in the literature [32,[54][55][56][57][58][59], leads to a static inconsistency at the boundary. The aim of the present work is to show that all higher-order theories written in the form of Eq. (3a) lead to an erroneous vanishing of the shear force at clamped edges and an underprediction of transverse deflection for large thickness to side ratios and orthotropy ratios. Previously, Levinson [20] touched upon this behaviour by stating that at a clamped end the average rotation of the cross-section can not be specified, while Tessler et al. [52] discussed the issue for different Zig-Zag theories but did not elucidate the matter in detail.
As shown in Fig. 2 the analysis presented here pertains to a plate that is infinitely wide in the transverse y-direction, with finite length L in the x-direction and support conditions defined along the infinitely wide edges x ¼ x A and x ¼ x B . This configuration is chosen to reduce the structural complexity by enforcing all unknown fields to be functions of the x and z coordinates only, i.e. reducing the analysis to a beam in plane strain. As a result, the discussion only deals with the essential boundary condition dw ;x ¼ 0 at edges x ¼ x A and x ¼ x B allowing for simpler derivations and clearer discussion of the arguments made. However without loss of generality, all comments apply to two-dimensional plates and shells because the inclusion of the lateral displacement field u y in the principle of virtual displacements introduces the essential boundary condition dw ;y ¼ 0 at clamped edges y ¼ y A and y ¼ y B .
The rest of the paper is structured as follows. Section 2 defines the problem that is used to compare different theories. Section 3 compares the classical FSDT of Mindlin with an alternative formulation that includes the CPT rotation w ;i . In particular, the effects on modelling clamped boundary conditions in a physically and mathematically correct manner is discussed. Section 4 extends the analysis to higher-order theories and shows how Reddy's third-order theory results in a vanishing shear force at a clamped edge and associated underpredictions of transverse displacement u z compared to high-fidelity 3D Finite Element solutions. In Section 5 the concept of general higher-order theories along the lines of Carrera's Unified Formulation and the Generalized Unified Formulation is introduced, which are shown to overcome the deficiencies of higher-order theories that include the w ;i term. In fact, it is demonstrated that shear strains automatically vanish at the top and bottom surfaces when the order of the generalised theory is sufficient to capture all higher-order effects, obviating the need for enforcing this boundary condition a priori. A metric for assessing the order of the expansion required is developed. Finally, the conclusions of the present work and suggestions are drawn in Section 6.

Problem definition
For simplicity, but without lack of generality, the following analysis considers an infinitely wide, 0 orthotropic layer in cylindrical bending caused by an arbitrary transverse loading of magnitude q 0 (Fig. 2). The plate has thickness t and axial length L and may take any support condition at the two boundary edges x ¼ x A and x ¼ x B which are henceforth referred to as the ends of the plate. The plate may consequently be analysed using the plane strain condition and the structural behaviour reduced to a function of the axial coordinate x and transverse coordinate z. These assumptions are made such that factors driving the governing mechanics can be elucidated in a straightforward manner. For all problems solved, the deflection and stress results are stated as normalised quantities, w ¼ Therefore absolute magnitudes of material properties play no role; only the orthotropy ratio E 11 =G 13 , thickness to length ratio t=L, the Poisson's ratios v 12 ¼ v 13 ¼ 0:25 and the boundary conditions are of significance. A 3D FEM model in the commercial software package ABAQUS was used as a benchmark to validate the results. For all cases considered the plate was discretised with 400 and 200 linear C3D8R elements along the axial length and through the thickness, respectively. To enforce the plane strain condition in the lateral direction a single C3D8R element was applied along the width and lateral expansion prevented. The ensuing model of 80,000 elements yields converged results to within 0.1% for all cases presented in Sections 3-5. Applied external loading is split equally between the top and bottom surfaces of the plate to minimise local deformations.
All analytical formulations presented here are solved in the strong form using an implementation of the Differential Quadrature Method (DQM) in MATLAB. DQM is a numerical technique  in the 1970s by Bellmann et al. [60] for solving differential boundary value problems. Since then DQM has been shown to give robust and efficient solutions to many problems in structural mechanics [61].

First-order shear theory
To begin, we analyse the cylindrical bending of the infinitely wide plate using the First-order Shear Deformation Theory (FSDT). One variant of FSDT is Mindlin's plate theory in which the functional unknowns are the mid-plane displacement u 0 , the transverse deflection w and the average rotation of the cross-section h. The average rotation of the cross-section may be assumed to comprise of the Kirchhoff bending rotation Àw ;x and the average shear rotation c. Thus, an alternative way of writing the in-plane displacement assumption of Mindlin's theory would be to replace h with Àw ;x þ c. In the following two examples, we investigate the effect this has on the governing field equations and boundary conditions as derived in a variationally consistent manner from the principle of virtual displacements (PVD).

Mindlin plate
For cylindrical bending, Mindlin's plate theory assumes axial and transverse displacements in the following form Using the kinematic relations between strains and displacements the axial strain x and transverse shear strain c xz are given by The principle of virtual displacements states that a body is in equilibrium if the virtual work done by the equilibrium forces, when the body is perturbed by a virtual amount dũ from the true configurationũ is zero. With regard to a plate subjected to cylindrical bending in the plane strain condition, the virtual work done by the virtual displacement dũ is dP ¼ where q is the net transverse pressure and S 2 is the boundary surface on which the stressesr x andŝ xz are prescribed. Substituting the strains of Eq. (8) into the PVD statement Eq. (9) yields dP ¼ The stresses r x and s xz are integrated through the thickness to define the stress resultants N; M and Q, where N is the in-plane load per unit width, M is the bending moment per unit width, Q is the shear force per unit width and k is the pertinent shear correction factor. The shear correction factor is needed to energetically account for the actual parabolic shear stress profile and is assumed to be equal to 5=6 [12]. Thus, Eq. (10) reduces to dP ¼ where x A and x B are the two supported ends of the plate. The governing equations and boundary conditions are derived using the calculus of variations. Integrating by parts those variational terms that feature derivatives we arrive at The Euler-Lagrange equations that cause the integral expression and the expression evaluated at x A and x B to vanish give the governing field and boundary equations, respectively. These governing field equations are while the essential and natural boundary conditions are Modelling a clamped boundary condition in Mindlin's plate theory is relatively straight forward. If we assume both ends x A and x B are rigidly built-in there cannot be any in-plane or transverse displacement and the plate cross-section can not rotate. Thus, Note, this boundary condition does not specify that the rotational component w ;x ¼ 0 at the ends x A and x B . In fact we may have a non-zero value of w ;x due to the presence of transverse shear rotation. Furthermore, we know from basic equilibrium that a non-zero shear force Q is required at the clamped edges. Using the definition of Q in Eq. (11) with the transverse shear constitutive equation and kinematics we have Thus, for a clamped edge with h ¼ 0 the value of w ;x defines the magnitude of the shear force at the support. In Kirchhoff's plate theory the shear rotation is assumed to be zero with the effect that clamped boundary conditions need to be enforced by setting w ;x ¼ 0. Physically, this boundary condition is reached asymptotically as the thickness-to-length ratio t=L of the plate approaches zero. Because all plates have finite thickness and thus finite shear deformation, imposing the condition w ;x ¼ 0 is physically incorrect and leads to large inaccuracies as the thickness increases, as shown later.

Alternative first-order shear theory
The observations of the previous section may now be compared to a theory where the average rotation of the cross-section h is replaced with a sum of the Kirchhoff rotation w ;x and the shear rotation c. Under these circumstances the displacement field is given by which leads to a new strain field, Repeating the same variational analysis of the previous section using the PVD definition Eq. (9) and the new strain field Eq. (19) results in where N; M and Q are as previously defined in Eq. (11). The governing field equations derived from the integral expressions are while the essential and natural boundary conditions at x A and x B are The same redundancy is also shown in the boundary conditions. The natural boundary conditions in Eqs. (22b) and (22c) are the same. Solving the boundary value problem of differential equilibrium Eqs. (21) in a mathematically consistent manner requires four boundary conditions at each end, i.e. eight boundary conditions in total. For simply-supported boundary conditions where M ¼ 0 at either end we thus only have six boundary conditions to apply due to this redundancy. It is of course still possible to solve the governing equations by assuming mode shapes for u 0 ; w and c that satisfy the boundary conditions. However, the problem in itself is not mathematically well-determined.
In fact, if we consider clamped boundary conditions we fail to solve the problem at all. If both x A and x B are both rigidly built-in the boundary conditions u 0 ¼ w ¼ 0 need to be satisfied. The bending moment M is non-zero at both ends such that a value for both c and w ;x needs to be prescribed on the boundary. For a generic load condition q the values of both c and w ;x are unknown as both c and w are variables of the theory. Thus, the only condition we can apply to adequately constrain the boundary value problem is c ¼ w ;x ¼ 0.
This statement causes two anomalies. First, we know that w ;x is in fact non-zero at the boundary due to the presence of transverse shear rotation. Second, if c ¼ 0 at x A and x B then the shear force vanishes at the boundary, where the shear strain c xz ¼ c according to Eq. (19b). On the contrary, from simple equilibrium it is clear that the shear force is finite at the supports x A and x B .

Higher-order theories featuring Kirchhoff rotation
In FSDT the presence of the Kirchhoff rotation w ;x has introduced two critical inconsistencies. The argument presented in the previous section is now extended to higher-order theories. Reddy's third-order theory is used as an example but the observations apply to any higher-order theory that is written in the general form of Eq. (3a). These include, but are not limited to, the theories of Ambartsumyan [22], Touratier [27], Soldatos [30] and Karama [32].

Reddy's higher-order theory
To derive Reddy's third-order theory we start with a cubic expansion of the in-plane displacement field, where h is the average rotation of the cross-section and f and n are higher-order rotations. The tractions at the top and bottom surfaces z ¼ AEt=2 are known a priori and this boundary condition is enforced on the in-plane displacement field. The transverse shear from kinematics is given by Assuming zero shear traction at the top and bottom surfaces z ¼ AEt=2 the transverse shear strain must vanish accordingly, Thus, the second-order rotation f is eliminated from the inplane displacement expansion Eq. (24a) due to symmetry, and the third-order rotation n is replaced by a term involving the average rotation h and Kirchhoff rotation w ;x . The modified displacement field reads and the strain field is given by The governing field and boundary equations are derived using the principle of virtual displacements Eq. (9) and the new strain field Eq. (28), where c 1 ¼ 4 3t 2 and c 2 ¼ 4 t 2 , and N and M are defined as previously in Eq. (11). The higher-order moment P, shear force Q and higherorder shear force R are given by Note that the shear correction factor k is no longer required in the definition of the shear force Q because a parabolic shear strain profile has been assumed. The governing field and boundary conditions are the Euler-Lagrange equations of the virtual displacement statement Eq. (29). The governing field equations from the integral expressions are and the boundary conditions at the ends x A and x B are First, it is noticed that the force boundary condition on the higher-order moment P in Eq. (32c) also features in Eq. (32b). Thus, we face the same problem with uniqueness of boundary conditions described in the previous section that is required to have a welldefined boundary value problem. For simply supported boundary conditions, and it is possible to define mode shapes that satisfy the boundary conditions exactly and then solve for the unknown coefficients. However, when solving a problem with built-in supports the situation is more difficult. The conditions are readily applied at either end as there can be no in-plane movement u 0 , no transverse deflection w and no average rotation h of the cross-section. This leaves the third boundary condition Eq. (32c) where either the higher-order moment P or the Kirchhoff rotation w ;x need to be defined. Both the bending moment M and higher-order moment P are non-zero at the supports because of a unique axial stressr x at this location. Furthermore, for a generic distributed transverse load q the value ofP at the supports is not known a priori. Therefore, we are forced to define a displacement boundary condition on the Kirchhoff rotation w ;x . Reddy applies the following boundary conditions for a clamped edge,  [32,54,55,[57][58][59] for various other higherorder theories as well. As described in Section 3, the boundary condition w ;x is physically inaccurate due to the presence of shear rotation at a clamped edge. Consider, for example, the transverse deflection plots of a plate with thickness-to-length ratio t=L ¼ 1 : 8 and orthotropy ratio E 11 =G 13 ¼ 50 clamped at the two ends and loaded with a uniformly distributed pressure. Fig. 3 shows the discrepancy of the rotation at the support when compared to Mindlin's theory and 3D FEM. While w ;x is non-zero in Mindlin's theory and in 3D FEM due to the presence of transverse shear, the rotation w ;x is forced to vanish in Reddy's theory. As a result, the plot shows that the boundary condition w ;x ¼ 0 overconstrains the structure and leads to stiffer behaviour. Mindlin's theory overpredicts the 3D FEM solution because for the chosen configuration higher-order effects are important. Mindlin's theory fails to capture the ''stress-channelling'' of axial stress r x towards the surfaces, thereby underpredicting the maximum stress and overpredicting the strain, i.e. the transverse deflection.
Furthermore, enforcing boundary condition Eq. (35) causes the shear force Q and higher-order shear force R to erroneously vanish at the supports, Thus, Reddy's higher-order theory leads to an inconsistency at a clamped edge between the shear forces obtained from constitutive and equilibrium equations. To overcome this inconsistency a nonzero value of w ;x could be defined. Alternatively, a boundary condition on w ;x could be not defined at all. The former is not possible as the slope w ;x depends on the loading condition, plate dimensions and material properties and is thus a quantity to be determined as a result of solving the problem. The latter is infeasible as an additional boundary condition is required to properly define the boundary value problem.

Shear force in Reddy's higher-order theory
Consider a cantilevered plate of thickness-to-length ratio t=L ¼ 1 : 10 and orthotropy ratio E 11 =G 13 ¼ 25 subject to a transverse shear traction q 0 at the free end. From simple equilibrium of forces and moments we know that the shear force must be constant along the length of the plate. The variation of the shear force from both constitutive and equilibrium conditions along the span of the plate, as derived from Reddy's plate theory, is shown in Fig. 4(a). It is apparent that at the built-in support x ¼ 0 there is a discrepancy between the constitutive shear force Q of Eq. (30) and the shear force V ¼ M ;x derived from equilibrium. Furthermore, the constitutive shear force can be seen to converge from zero to the correct value of unity some distance away from the built-in support. Fig. 4(b) shows an equivalent plot for a plate that is clamped at both ends and loaded by a uniformly distributed pressure of magnitude q 0 . In this case the constitutive shear force vanishes at both supports and converges to the linearly varying distribution derived from equilibrium some distance away from the supports. In both cases the convergence distance depends on the magnitude of the thickness-to-length ratio t=L and the orthotropy ratio E 11 =G 13 , as these ratios are instrumental in quantifying transverse shear flexibility. The variation of the constitutive shear force along the length of the cantilevered plate for different thickness-tolength ratios t=L is shown in Fig. 5(a) and for different orthotropy ratios E 11 =G 13 is shown in Fig. 5(b). As t=L and E 11 =G 13 increase so does the distance required to converge. For t=L ! 0 and E 11 =G 13 ! 0 the convergence distance approaches zero because the plate approaches the idealised condition of pure bending with w ;x ¼ 0. However, this trend is asymptotic such that the shear force at the support condition is always equal to zero for any finite thickness-to-length ratio t=L and orthotropy ratio E 11 =G 13 .
Also note that Mindlin's plate theory does not lead to this inconsistency and thus results in the same shear force along the length of the plate for both constitutive and equilibrium derivations. This is because in Mindlin's plate theory no constraint is placed on w ;x at the boundary.
The physically incorrect boundary condition w ;x ¼ 0 does not only lead to an anomaly for the shear force but influences the bending deflection and stress results as a whole. Because the plate may in fact shear at the boundary, the boundary condition w ;x ¼ 0 overconstrains the structure and leads to unconservative results for transverse deflection. The inaccuracy increases with increasing thickness-to-length ratio t=L and increasing orthotropy ratio E 11 =G 13 . A comparison of normalised transverse deflection w at the midspan for both Mindlin (FSDT) and Reddy's higher-order plate theory (RHOT) against the high fidelity 3D FEM solution for different thickness-to-length ratios t=L is shown in Table 2. The same comparison for varying orthotropy ratio E 11 =G 13 is shown in Table 3. The normalised solution of 3D FEM is stated explicitly and the errors in FSDT and RHOT are given as percentages. The solutions are calculated for a plate clamped at both ends and loaded by a uniformly distributed pressure q 0 across the span L. Table 2 shows that the accuracy of w for FSDT and RHOT are similar for small values of t=L. However, as the thickness increases Reddy's higher-order theory is less accurate than Mindlin's plate theory. The results show that RHOT always leads to an underprediction of w compared to 3D FEM which arises due to the stiffening effect of w ;x ¼ 0 at the boundary. For FSDT the error is always on the conservative side. Furthermore, for t=L ¼ 1 : 5 the 9.44% error of RHOT is more than three times the magnitude of the FSDT error. Finally, the results for w in Table 3 show a similar trend of increasing inaccuracy for both FSDT and RHOT as the orthotropy ratio E 11 =G 13 increases, with the error of RHOT generally two times greater than that of FSDT.
The errors in FSDT may largely be attributed to a failure in capturing higher-order effects that become important as t=L and E 11 =G 13 increase. RHOT on the other hand captures higher-order effects due to the cubic displacement formulation and the errors stem largely from the erroneous boundary condition w ;x ¼ 0.

General higher-order theories
Based on the previous findings a formulation is sought that: 1. Captures higher-order effects. 2. Allows transverse shear stresses to vanish at the surfaces. 3. Gives a meaningful, non-zero shear force at a clamped edge, i.e.
Kirchhoff's term w ;x does not appear in the axial displacement expansion.
To achieve this a general higher-order theory following Carrera et al. [9,42,43] is presented.

Model derivation
The most expedient way to achieve the above conditions is by a natural extension to Mindlin's theory; that is, the axial displacement field u x is expanded in generalised form. For example u x may be written as a linear combination of a power series in z and an unknown displacement field U, where h is the average rotation of the cross-section, and f; n; w are parabolic, quartic and sixth-order distortions of the cross-section, respectively. Eq. (37) is based on the same idea as the generalised expansion in Carrera's Unified Formulation and the Generalized Unified Formulation, where the order of a theory is defined by the highest order of the transverse coordinate z. Note that even terms in z are unnecessary due to the intrinsic symmetry of the layer considered.
The generalised strain fields are given by The shear strain is not forced to vanish at the top and bottom surfaces a priori, such that the shear strain and stress profiles are of higher-order but are free to ''float'' away from zero at the surfaces. In the following, it is shown that the shear strain automatically vanishes at the surfaces if the order of the expansion adequately captures the structural behaviour. The greater the value of the thickness-to-length ratio t=L and orthotropy ratio E 11 =G 13 the more terms in the expansion are required to achieve this. Thus, enforcing the shear strain to vanish a priori, as in RHOT, is not necessary and the inconsistency due to the w ;x term can be prevented.
Using the principle of virtual displacements Eq. (9) with new strain field Eq. (38) gives dP ¼ Integrating Eq. (39) in the z-direction by defining the following stress resultants where T signifies the transpose of the vectors N and Q, and then integrating by parts gives Note, that no shear correction factor is required for Q and Q due to the presence of higher-order terms. The set of governing field equations is derived by setting the integral expression to zero, while the essential and natural boundary conditions at x A and x B are   Table 3 Normalised transverse displacement w of 3D FEM results compared to Mindlin's (FSDT), Reddy's higher-order (RHOT), third-order (3HOT) and fifth-order (5HOT) solutions, for a plate clamped at two ends and loaded by a uniformly distributed pressure.
Normalised transverse deflection w at x=L ¼ 0:5 Thus, we have arrived at a generalised set of governing field equations and boundary conditions. The number of equations in the set corresponding to dU depends on the order of the axial displacement expansion, while the transverse equilibrium equation of dw is always fixed. At a clamped boundary condition we can set all rotations and the transverse deflection to be zero, i.e. U ¼ w ¼ 0. This formulation allows the support shear force to be finite at a clamped edge. Also note that the shear stress at a clamped edge s xz ¼ G xz w ;x is independent of the higher-order field F ;z and is therefore a constant. This result agrees with the physical reality at a built-in support, whereby the whole cross-section is sheared by the same amount. Thus the shear stressŝ xz at the built-in support is constant and non-zero throughout the whole thickness. Tables 2 and 3, previously used to compare the results of normalised transverse deformation w at the midspan for FSDT and RHOT, also include the results of cubic (3HOT) and quintic (5HOT) generalised theories. As expected, the fifth-order theory provides slightly more accurate results than the third-order theory for both increasing thickness-to-length ratio t=L and orthotropy ratio E 11 =G 13 . Also, both 3HOT and 5HOT considerably improve on the accuracy of Reddy's higher-order theory, as the inconsistency due to the w ;x boundary condition is removed. The difference is especially striking for E 11 =G 13 ¼ 200 in Table 3 where the error of 3HOT and 5HOT are one and two orders of magnitude smaller than RHOT, respectively.

Comparison with FSDT and RHOT
The maximum normalised axial stress at the midspan of the plate r x as predicted by FSDT, RHOT, 3HOT and 5HOT are compared against 3D FEM in Table 4 for increasing thickness-to-length ratio t=L. FSDT significantly underpredicts the maximum stress for large t=L values. This behaviour is to be expected as the linear stress assumption can not capture the higher-order ''stresschannelling'' effect for large t=L ratios. Even though the boundary condition on w ;x in RHOT leads to inaccurate transverse deflection results, the stress solutions maintain good accuracy for all t=L presented. Thus, it seems that strain and stress results, being based on derivatives of the displacements, are not affected to the same degree as the displacements themselves. Nevertheless, both 3HOT and 5HOT always outperform Reddy's higher-order theory.
A similar trend is shown in Table 5 for increasing orthotropy ratio E 11 =G 13 . As the orthotropy ratio increases the plate becomes relatively more flexible in shear which increases the distortion of the plate cross-section. This effect increases higher-order ''stress-channelling'' effects. Thus, the stress profile through the thickness transitions from pre-dominantly linear, to cubic, to quintic and to further higher-order fields. Both 3HOT and 5HOT considerably improve upon FSDT and RHOT, giving nominal errors. For increasing E 11 =G 13 the axial stress r x from RHOT now show errors up to 8%. This amount is more than four times the error in 3HOT and 5HOT.

Hierarchical Modelling
One of the advantages of RHOT is that the transverse shear stresses are forced to vanish at the top and bottom surfaces. This is a physically meaningful result that at the same time reduces the number of variables to that of FSDT. However, enforcing this constraint introduces the Kirchhoff rotation w ;x into the formulation and leads, as we have seen, to underpredictions in transverse displacement and vanishing of the support shear force for clamped boundary conditions. In 3HOT and 5HOT the transverse shear strain is not enforced to vanish on the top and bottom surfaces. Nevertheless, if the chosen order of the theory adequately captures the structural behaviour, the transverse shear stress naturally disappears at the top and bottom surfaces.
As an example, consider a plate of thickness-to-length ratio t=L ¼ 1 : 10 and orthotropy ratio E 11 =G 13 ¼ 50. The plate is simply supported at either end and loaded by a sinusoidally distributed pressure q ¼ q 0 sin px L À Á . This loading configuration is chosen as it allows the through-thickness results to be compared to Pagano's 3D elasticity solution [7].
In Fig. 6(a) the normalised axial stress r x at the midspan of the plate is plotted through the plate cross-section for Pagano, FSDT, 3HOT and 5HOT solutions. The linear profile of FSDT fails to capture the higher-order field but both the 3HOT and 5HOT are accurate throughout the whole plate cross-section. As a result, the normalised transverse shear stress profile s xz at the support x ¼ 0 is accurately captured by both 3HOT and 5HOT ( Fig. 6(b)). For 5HOT the shear stress vanishes exactly at the top and bottom surfaces while for 3HOT a small residual remains. Now consider a plate t=L ¼ 1 : 8 and E 11 =G 13 ¼ 200 under the same loading conditions. Fig. 7(a) shows that the ''stresschannelling'' of r x towards the outside surfaces is more pronounced.
While 5HOT remains close to Pagano's solutions throughout the whole thickness there are some inaccuracies for 3HOT. This is because fifth-order effects are significant for this plate configuration. As a result, there is a considerable residual in s xz at the surfaces for 3HOT, while the solution for 5HOT remains accurate through the thickness and vanishes on the surfaces (Fig. 7(b)).
The authors believe that this behaviour is a direct result of the minimisation technique employed in the principle of virtual displacements. In essence, each theory captures the total potential energy through the volume of the body in an average sense (volume integral) as permitted by the order of the theory. If the true 3D behaviour of the structure is governed by higher-order behaviour, as in Fig. 7, the linear profile of FSDT must underpredict in Table 4 Normalised axial stress rx of 3D FEM results compared to Mindlin's (FSDT), Reddy's higher-order (RHOT), third-order (3HOT) and fifth-order (5HOT) solutions, for a plate clamped at two ends and loaded by a uniformly distributed pressure.   some parts and overpredict in others in order to arrive at a similar magnitude of total potential energy. The same argument holds for the transverse shear stress of 3HOT in Fig. 7(b). The quadratic shear stress expansion can not capture the actual higher-order behaviour accurately. Thus, it overpredicts at the surfaces and the midplane, while underpredicting at other points throughout the crosssection. This effect is essentially removed in 5HOT by allowing enough degrees of freedom to accurately capture the higher-order stress profile. The adequacy of a higher-order theory in modelling the structural behaviour can be ascertained by analysing the residual of the transverse shear stress at the surfaces. The energy associated with this residual is given by If this residual is, say, three orders of magnitude smaller than the average transverse shear energy through the thickness, then the error associated with this residual can be assumed to be negligible. Consequently, if Àt=2 s xz c xz dz   in Table 6. The result in Table 6 supports the qualitative observations made regarding Fig. 7(b) that 3HOT is inadequate in capturing the higher-order effects in the plate with t=L ¼ 1 : 8 and E 11 =G 13 ¼ 200. This result has been underlined in the table.

Asymptotic Expansion
It was shown in the previous sections that the structural behaviour is a function of both the orthotropy ratio E 11 =G 13 and the thickness to length ratio t=L. It is possible to combine these two factors into a single metric which governs the structural behaviour of the plate.
Consider an infinitely wide plate as depicted in Fig. 2 simply supported at either end and loaded by a sinusoidally distributed pressure q ¼ q 0 sin px L À Á . The governing Eqs. (14b) and (14c) for FSDT written in terms of w and h are dh : Dh ;xx À kGðw ;x þ hÞ ¼ 0; ð47aÞ where k is a pertinent shear correction factor, D is the bending rigidity and G the shear rigidity of the single layer defined by The ad-hoc assumptions satisfy the boundary conditions exactly. Substituting Eq. (49) into the governing Eqs. (47) and solving for the unknown coefficients w 0 and h 0 gives Substituting Eq. (48) into Eq. (50) yields where k ¼ Q 11 t L À Á 2 is a lay-up dependent ratio that governs the influence of transverse shear deformation. The origin of this parameter is discussed by Hu et al. [62] and has been used extensively in the literature to assess the effect of transverse shear deformation on the structural behaviour of beams, plates and shells. As FSDT only captures first-order effects, the order of k in Eq. (51) is consequently unity. A similar analysis is conducted for 3HOT. Thus, writing the governing Eqs. (42) in terms of the unknown variables h; f and w gives dh : Dh ;xx þ Ef ;xx À Gðw ;x þ hÞ À Hf ¼ 0; ð52aÞ where D and G are as previously defined in Eq. (48), F and I are the higher-order bending and transverse shear rigidities, respectively, and E and H are the first-order/higher-order coupling bending and transverse shear rigidities, respectively. These are defined by The assumptions satisfy the boundary conditions exactly. Substituting Eq. (54) into the governing Eqs. (52) and solving for the unknown coefficients w 0 ; h 0 and f 0 gives Substituting Eq. (53) into Eq. (55) yields which shows that the transverse displacement is now a higherorder function of k. The expression in Eq. (56) may be expanded using the binomial expansion and compared to the FSDT expression in Eq. (51). Thus, expanding Eq. (56) as a power series gives which shows that if the shear correction factor for FSDT is chosen to be k ¼ 5=6, i.e. the original value found by Reissner [45], then the first-order coefficients of k in Eqs. (51) and (57) are equal. Thus, this shear correction factor of k ¼ 5=6 only attempts to impose the firstorder structural effects associated with k to that of FSDT. The shear correction factor may be generalised to include all higher-order terms of k found in 3HOT by equating Eqs. (51) and (56) 1 þ The expression in Eq. (58) shows that the shear correction factor k ! 5=6 as k ! 0. Furthermore, as k ! 1 the value of k ! 1. This suggests that the shear correction factor k ¼ 5=6 pertains to a fictitious layer of infinitesimal thickness or infinite length with perfectly linear and parabolic variations of r x and s xz through the thickness, respectively. As k increases and the ''stress-channelling'' effect becomes more significant the stresses r x and s xz transition to increasingly higher-order profiles, as previously shown in Figs. 6 and 7. As k increases and the in-plane stress r x is channelled towards the outside surfaces the transverse shear stress s xz is more evenly distributed at the centre, as shown in Fig. 8; an effect which is analogous to the behaviour observed in a sandwich beam with stiff face layers and a compliant core. With more of the crosssection sheared by the same amount, the shear correction factor increases because the energetic difference between the constant shear stress profile of FSDT and the actual profile is decreasing. Thus, in the limiting case of constant shear stress through the thickness the shear correction factor approaches unity. In essence, Eq. (58) generalises the shear correction factor for a plate of finite thickness to any order of k to the accuracy provided by 3HOT. The same process may readily be extended to other higher-order theories in order to correct FSDT to capture global effects predicted by 5HOT and 7HOT. The corresponding magnitudes of transverse deflection w 5HOT 0 and w 7HOT 0 are shown in Appendix A.
In Table 7 the normalised transverse deflection from FSDT calculated using the generic shear correction factor and the classical value of k ¼ 5=6 are compared. The table shows that the percentage error with respect to Pagano's 3D elasticity solution [7] is significantly reduced when the shear correction factor in Eq. (58) is used.
The local axial stress results of the two FSDT theories however, are unchanged because the through-thickness assumption remains linear such that higher-order ''stress-channelling'' effects cannot be captured. Using the solutions for the displacement variables of FSDT in Eq. (50) and the kinematic and constitutive relations the axial stress field for FSDT is given by which is unchanged from the solution of CLA. Similarly using the solution for the displacement variables of 3HOT in Eq. (55) we get Expanding Eq. (60) as a power series, shows that the higher-order solution modifies the FSDT axial stress field in two ways. First, the average slope of the through-thickness field is modified by the power series coefficient in the first bracket of Eq. (61). Second, a higher-order component that governs the extent of ''stress-channelling'' is given by the power series coefficient of ðz=tÞ 2 . In this particular case the higher-order component has a leading coefficient that is 20=3 times greater than the leading coefficient of the average slope change. Thus, including the cubic term in the displacement field contributes more to the axial stress field through the introduction of the higher-order ''stress-channelling'' effect than by correcting the linear component.

Conclusions
Static inconsistencies in modelling clamped boundary conditions using displacement-based, axiomatic, higher-order theories have been discussed. Enforcing the boundary condition of vanishing transverse shear strain at the top and bottom surfaces a priori introduces the Kirchhoff rotations w ;x and w ;y into the expansion for u x and u y , respectively. If the governing differential equations are derived using the principle of virtual displacements an essential boundary condition on w ;i perpendicular to an edge arises which is forced to vanish to properly constrain the boundary value problem at a clamped edge. This condition is physically inaccurate as the plate can actually rotate at a clamped edge due to the presence of transverse shear rotation. Furthermore, this condition on perpendicular rotation w ;i ¼ 0 causes the shear force, as derived from constitutive equations, to vanish at the clamped edge. Such a condition is erroneous when compared to simple transverse equilibrium conditions. Finally, preventing rotation at the clamped edge overconstrains the structure and leads to underpredictions of transverse deflections and overpredictions of axial stresses.
When the in-plane displacement fields are written as a general power series in terms of the transverse coordinate z, as proposed in Carrera's Unified Formulation, this inconsistency does not occur. Furthermore, if the order of the theory is sufficient to capture all higher-order effects the transverse shear stresses automatically vanish at the top and bottom surfaces, obviating the need to apply this constraint a priori. Based on this insight a non-dimensional parameter based on the transverse shear strain energy at the surfaces was introduced to gauge the accuracy of a higher-order theory. Finally, it was shown that the structural behaviour of a single layer plate in bending is a function of the parameter k ¼ Q 11 The parameter k can be used to derive shear correction factors for FSDT that allow the transverse bending deflection results to match any higher-order theory.    1) shows that the higher-order theories progressively correct higher-order terms in k. As a result the expansion of 3HOT converges to k 2 , 5HOT to k 4 and 7HOT to k 6 etc. compared to the infinite solution. Even though the expression for w 0 can be seen to follow a decaying progression, no closed form holonomic sequence could be determined. Table A.8 shows the convergence of normalised transverse deflection w 0 Dp 4 q 0 L 4 with increasing order of k n up to n ¼ 5 for three different materials with large thickness-to-length ratio t=L ¼ 0:3. The three materials are a metalic isotropic material with Poisson's ratio t ¼ 0:3, an industrial grade carbon-fibre pre-preg like IM7 8552 and a highly orthotropic lamina. The results show that for metalic isotropic materials the solution converges for n ¼ 1, i.e. a first-order theory solution. For materials with k % 3 such as a very thick IM7 8552 composite the solution converges for n ¼ 4, i.e. a third-order theory solution. For the highly orthotropic material the full series solution between 5HOT and 7HOT are not much different but a power series solution in terms of k n needs to be expanded past n ¼ 5 to achieve convergence.

Acknowledgments
Finally, the authors wish to point out that higher-order solutions past 5HOT may potentially be meaningless. The structural behaviour for large values of k where these higher-order terms become significant may exhibit greater proportion of transverse normal deformation, which has been ignored in the present analysis. In these situations a full 3D solution is required. q 0 L 4 with increasing order of k n for three different materials with thickness-to-length ratio t=L ¼ 0:3. Convergence is compared to the full series solution (n ! 1) of 3HOT, 5HOT and 7HOT.