Static-kinematic duality in beams, plates, shells and its central role in the finite element method

Abstract Static and kinematic matrix operator equations are revisited for one-, two-, and three-dimensional deformable bodies. In particular, the elastic problem is formulated in the details in the case of arches, cylinders, circular plates, thin domes, and, through an induction process, shells of revolution. It is emphasized how the static and kinematic matrix operators are one the adjoint of the other, and then demonstrated through the definition of stiffness matrix and the application of virtual work principle. From the matrix operator formulation it clearly emerges the identity of the usual Finite Element Method definition of elastic stiffness matrix and the classical definition of Ritz-Galerkin matrix.


Introduction
The static-kinematic duality leads to a simple and direct demonstration of the Principle of Virtual Work for deformable bodies, and vice-versa [2][3][4][5]. The two concepts imply each other. Such a demonstration derives from the representation of the elastic problem in a symmetrical manner by combining the three fundamental relationsindefinite equations of equilibrium, kinematic equations as definition of deformation characteristics, and constitutive equations -in a single matrix operator equation where the unknown is represented by the displacement vector.
The definition of the static and kinematic matrix operator equations for 1-, 2-, 3-D deformable bodies and the formulation of the elastic problem for arches, cylinders, circular plates, thin domes, and shells of revolution will herein be presented. It is possible to verify that each kinematic matrix operator is the adjoint of the corresponding static matrix operator, and vice-versa.
From the formulation of the Finite Element Method stiffness matrix [1,10], clearly emerges its identity with the Ritz-Galerkin matrix. After more than half a century from the pioneering proposal of the Finite Element Method, it appears to be the right time today to reconsider its foundations, to emphasize its strict connections with the more general Ritz-Galerkin Method, and, based on the statickinematic duality, to reformulate it in a more modern way. We can not exclude that such a new vision could pave the way to more advanced and simple discretization methods.

Three-dimensional body
It is considered an elementary parallelepiped with the sides parallel to the coordinate axes, of length dx, dy, dz, respectively. On the opposite faces of the parallelepiped there act components of stress which, but for an infinitesimal increment, are equal to one another. Equilibrium to translation in the X direction, for instance, imposes ∂σx ∂x dx(dydz) + ∂τyx ∂y dy(dxdz) + ∂τzx ∂z dz(dxdy) (1) where only the increments of stress, multiplied by the elementary areas on which they act, and the body force, multiplied by the elementary volume in which it acts, are present. Dividing Equation (1) by the elementary volume dV = dx dy dz, we obtain the first of the indefinite equations of equilibrium: The analogous equations of equilibrium in the Y and Z directions appear as follows: ∂τxy ∂x + ∂σy ∂y + ∂τzy ∂z + Fy = 0 (3) The previous equations may also be obtained by integration.
We are now able to express in matrix form the systems of differential equations that govern, on the one hand, equilibrium and, on the other, congruence. These two systems have intimately connected formal structures, as we shall show in this section.
As regards equilibrium, the indefinite equations of equilibrium may be reproposed, in matrix form, with the stress components ordered in a stress vector premultiplied by a matrix operator (3 × 6): The static equations can be written in compact form As regards congruence, the six independent components of deformation can be ordered in a strain vector, which can be obtained by premultiplying formally the displacement vector by a (6 × 3) matrix operator: The relations (7) can also be written in compact form, and are called kinematic equations. The purely differential static operator is exactly the transpose of the kinematic one appearing in Equation (7).
In this way, just as in the mechanics of rigid bodies, the static matrix is the transpose of the kinematic one [8], also in the mechanics of deformable bodies there exists the same profound interconnection between the two matrix operators. In the case of the mechanics of rigid bodies, we know how this interconnection implies the validity of the Principle of Virtual Work. The validity of this principle can, on the other hand, be extended to the case of deformable bodies, precisely on the basis of the static-kinematic duality previously emphasized.
The three indefinite equations of equilibrium do not suffice to determine the six components of stress. On the other hand, by adding to them the six elastic constitutive equations, we obtain a system of nine differential equations in nine unknowns: σx, σy, σz, τxy, τxz, τyz; u, υ, w.
In matrix form, it is possible to give a very synthetic and expressive representation of the linear elastic problem, considering as the primary unknown the displacement vector {η}. If in the static equation we introduce the constitutive law, and then the kinematic equation, we obtain a matrix operator equation, called Lamé equation: The matrix and second order differential operator in round brackets is called the Lamé operator, It turns out to be a (3 × 3) matrix and, in non-homogeneous problems, where the Hessian matrix [H] is a function of the point, it too is a function of the point.
Recalling the boundary equations of equivalence and assuming that they hold good on a portion Sp of the external surface of the body and that, on the complementary portion Sη, there is imposed a congruent field of displacements {η 0 }, the three-dimensional elastic problem can be synthesized as follows:

Beam with rectilinear axis
It is considered an elementary portion of a beam with rectilinear axis and a cross section that is symmetrical with respect to the Y axis. Let this portion be subjected to a bending moment Mx and to a shearing force Ty. Deformations due to these two characteristics will produce relative displacements between the centroids of the two extreme cross sections of the beam portion, exclusively in the transverse direction of the Y axis. In the case of the shear Ty, we have ( Fig. 1(a)): where dυ T is the relative displacement in the Y direction due to the shear, y is the shearing strain dual of the shearing force and dz is the length of the infinitesimal element of beam. In the case of bending moment, and considering the rigid rotation φx of the element, we have ( Fig. 1(b)): having neglected the infinitesimals of a higher order due to the curvature, i.e. to the slope variation dφx ( Fig. 1(c)): Summing up the two significant contributions of the shearing strain and of the rigid rotation, we obtain At this point we are able to formulate the fundamental equations of the elastic problem for beams with a rectilinear axis. The kinematic equations constitute, as in the case of the three-dimensional body, the definition of the characteristics of deformation as functions of the generalized displacements [7]: where among the components of the deformation vector appear the shearing strain y , the axial dilatation εz, and the flexural curvature χx, whilst among the components of the displacement vector appears, in addition to the ordinary components, υ and w, also the rotation φx. The transformation matrix is differential and shows on the diagonal the total derivative d/dz, whereas the off-diagonal terms are all zero except for one, which is equal to +1. The same relation may be written in compact form, where {q} indicates the vector of the deformation characteristics, {η} the vector of generalized displacements, and [∂] the kinematic matrix operator.
On the other hand, also the indefinite equations of equilibrium can be presented in matrix form [7]: where among the components of the vector of static characteristics appear the shearing force Ty, the axial force N, and the bending moment Mx, whilst among the components of the vector of external generalized forces appears, in addition to the transverse distributed load q and the axial distributed load, p, also the bending distributed moment mx. The matrix operator presents the total derivative d/dz in all diagonal positions, whereas the off-diagonal terms are all zero, except for one, which is equal to −1 and expresses the identity of the shear with the derivative of the bending moment (neglecting the distributed moment mx). The static matrix operator (21) is equal to the transpose of the kinematic matrix operator (19), but for the algebraical terms which present opposite signs. This is said to be the adjoint operator of the previous one, and vice versa. The reason why the algebraical terms change sign will emerge clearly afterwards.
In compact form we can write where {Q} is the vector of static characteristics, and {F} is the vector of external forces. It should be noted that, unlike in the case of the threedimensional body, the matrices [∂] and [∂] * are square (3 × 3). This means that it is possible to determine the internal reactions of a statically determinate beam by using only the equilibrium equations. This, unfortunately, does not occur in the case of a three-dimensional body, also if constrained isostatically, the stress field of which may be determined only by applying, in addition to the static equations, also the constitutive and the kinematic equations.
Having now at our disposal the equations of kinematics and statics and the constitutive equations for the beam, we can obtain Lamé equation in operator form: The matrix and differential operator of the second order in round brackets can be called the Lamé operator It turns out to be a (3 × 3) matrix and, in non-homogeneous problems, in which the Hessian matrix [H] is a function of the axial coordinate z, is also a function of z.
Finally, the boundary conditions may be conditions of equivalence at the ends, T is the identity matrix [1], with a value of unity corresponding to each differential term of the matrix [∂] * . The boundary conditions can, on the other hand, also represent displacements imposed at the ends: In conclusion, the elastic problem of the rectilinear beam can be summarized as follows: where the static boundary condition (28) holds good in the end point or points that are subjected to loading, as does the kinematic boundary condition (29) in the constrained end point or points. In the case where there are no conditions on the displacements, the loads {Q 0 } must constitute a self-balancing system.

Beam with curved axis
If we consider the element of beam with curved axis, the curvilinear coordinate s is considered as increasing as we proceed from left to right along the beam, while the angle dϑ is considered to be positive if it is counterclockwise. In accordance with the above conventions, also the radius of curvature r acquires an algebraic sign on the basis of the relation ds = rdϑ (30) Figure 2: Curved beam element: (a) reference system, (b) rigid rotation to be subtracted from the total one and due to a uniform tangential displacement, (c) additional axial dilation due to a uniform radial displacement.
As regards the generalized displacements of the generic cross section, the radial displacement υ is positive if it is in the positive direction of the Y* axis (where Y*Z* is a system of right-handed axes travelling along the axis of the beam), the tangential displacement w is positive if it is in the positive direction of the curvilinear coordinate s, and finally the variation of the angle φ is positive if it is counterclockwise ( Fig. 2(a)).
We shall now show how the kinematic equation must be modified to take into account the intrinsic curvature of the beam. The tangential displacement w produces, in fact, an apparent rotation φ(w), which is to be subtracted from the total rotation ( Fig. 2(b)): On the other hand, the radial displacement υ produces an additional axial dilation ε(υ) which is given by (Fig. 2(c)) As a consequence of an infinitesimal relative rotation dφ of the extreme cross sections of the beam element, the angle between the sections can be obtained as the sum (dϑ + dφ) of the initial and intrinsic relative rotation with the elastic and flexural relative rotation. The new curvature is then so that the variation of curvature is Substituting the relations (31), (32) and (34), which furnish respectively the rotation to be deducted, the additional axial dilation, and the variation in curvature, the kinematic equations for the curved beam appear as follows: The indefinite equations of equilibrium, or static equations, may be proposed in matrix form: It should be noted that, but for the algebraic signs of the non-differential terms, the static matrix is the transpose of the kinematic one. The rotation matrix which transforms the global reference system YZ into the local reference system Y*Z* is the following: so that the vectors of the external forces and of the generalized displacements in the local reference system may be expressed by premultiplying the respective vectors evaluated in the global reference system by the matrix [N] The static and kinematic equations can thus be expressed as follows: Substituting equation (43)  which is Lamé equation for curved beams, arches, and rings. The elastic problem can then be summarized as follows:

Cylindrical shell
Let us consider a cylindrical shell. The longitudinal coordinate along the generatrix is indicated with x, as well as the finite principal radius of curvature coincides with the radius R of the circular directrix (Fig. 3).
The kinematic equations are: where εx and ε ϑ are the longitudinal and circumferential dilations, respectively, x is the shearing strain, χx and χ ϑ are the longitudinal and circumferential curvatures, respectively, as well as u is the longitudinal displacement, w the transverse or normal displacement, and φx the rotation about the circumferences. On the other hand, the static equations are: where Nx and N ϑ are the longitudinal and circumferential unit forces, respectively, Tx is the longitudinal shearing force, Mx and M ϑ are the longitudinal and circumferential bending moments, respectively. We observe that the only external force acting on the shell is assumed to be a load q(x), normal to the middle surface. The variation in curvature χ ϑ vanishes, just as the moment M ϑ is not involved in any of the three equations of equilibrium.
The first of Equations (50) is the equation of equilibrium with regard to longitudinal translation, which gives Nx = constant. The second of Equations (50) is the equation of equilibrium with regard to radial or normal translation, while the third is the equation of equilibrium with regard to rotation about the parallel, Also in the case of the cylindrical shell the static and the kinematic matrix operators are the adjoint of one another.

Circular plate
The kinematic equations for a circular plate are: where εr and ε ϑ are the radial and circumferential dilations, respectively, r is the radial shearing strain, χr and χ ϑ are the radial and circumferential curvatures, respectively, as well as u is the radial displacement, w the transverse or normal displacement, and φx the rotation about the circumferences. On the other hand, the static equations are: where Nr and N ϑ are the radial and circumferential unit forces, respectively, Tr is the radial shearing force, Mr and M ϑ are the radial and circumferential bending moments, respectively.
Restricting the analysis to the flexural regime only and excluding the in-plane or membrane regime, we obtain the following equations, which are kinematic and static, respectively ( Fig. 4): The indefinite equations of equilibrium (57) represent a system of two differential equations in the three unknowns Tr, Mr, M ϑ . The polar symmetry thus reduces the degree of static indeterminacy of the deflected plate from two to one. The first of equations (57) represents the condition of equilibrium with regard to the normal translation of a plate element identified by two radii forming the angle dϑ, and by two circumferences of radius r and r + dr dTrrdϑ + Tr dr dϑ − qr dr dϑ = 0 (58) The second term of the foregoing equation is due to the greater length presented by the outermost arc of circumference. Dividing by the elementary area rdrdϑ, we obtain the first of equations (57). The second of equations (57) represents the condition of equilibrium with regard to rotation of the same plate element about the circumference of radius r: −Trrdϑdr + dMrrdϑ + Mrdrdϑ − M ϑ drdϑ = 0 (59) Also in this case the contribution of the third term is due to the greater length of the outermost arc of circumference.
In the case of the circular plate the matrix [∂] * can be obtained as before, by transposing and changing the sign of the algebraical terms of matrix [∂], but also by performing the following substitution [

Thin dome
Thin shells are shells of such small thickness that they present an altogether negligible flexural rigidity. These elements can sustain only compressive or tensile forces contained in the tangent plane.
As regards thin shells of revolution (thin domes), the kinematic and static equations simplify notably, since only forces along the meridians and the parallels, Ns and N ϑ , are present, as well as the displacements along the meridians and those normal to the middle surface, u and w, respectively (Fig. 5): The first of equations (61) represents the equilibrium to translation along the meridian, whereas the second represents the equilibrium to translation along the normal.
From the second of Equations (61) we obtain the fundamental algebraic relation which links the forces Ns and N ϑ while from the first we obtain the following differential equation: On the other hand, by means of equation (62) we can express N ϑ as a function of Ns, and this expression, inserted in equation (63), gives dNs ds + which is a differential equation with ordinary derivatives in the unknown function Ns(s

Shell of revolution
When a shell of revolution is loaded symmetrically with respect to the axis of simmetry Z, only the curvilinear coordinate s is present as an independent variable, while the displacement υ along the parallels vanishes, as well as the deformations sϑ , ϑ , χ sϑ , and the corresponding internal reactions N sϑ , T ϑ , M sϑ (Fig. 5): Observe that, again for reasons of symmetry, the conditions of equilibrium to translation along the parallel and to rotation about the meridian are identically satisfied and thus do not appear in Equation (67). Finally, we have three equations of equilibrium (respectively, with regard to translation along the meridian, to translation along the normal, and to rotation about the parallel) in the five static unknowns Ns, N ϑ , Ts, Ms, M ϑ . The elastic problem for shells of revolution thus has two degrees of internal redundancy, while the more general problem of shells with double curvature appears to have three degrees of internal redundancy.
Equations (67) are verified by imposing the above three conditions of equilibrium on an infinitesimal shell element, bounded by two meridians located at an infinitesimal angular distance rdϑ and by two parallels located at an infinitesimal distance ds.
The condition of equilibrium with regard to translation along the meridian yields the equation (Fig. 6(a) The condition of equilibrium with regard to translation along the normal n furnishes the equation (Fig. 6(c), 6(d), 6(e)): − Ns ds R 1 rdϑ − N ϑ ds dϑ cos α + dTsrdϑ + Tsdr dϑ (69) + q r ds dϑ = 0 which, divided by rdsdϑ, coincides with the second of equations (67). Finally, the condition of equilibrium with regard to rotation about the parallel furnishes the equation (Fig. 6(b), 6(c)): − Tsr dϑ ds + dMsr dϑ + Ms dr dϑ (70) − M ϑ sin α ds dϑ = 0 which, divided by rdsdϑ, coincides with the third of equations (67). Notice that, in the indefinite equations of equilibrium (67), five terms sinα/r are present. These contributions are due to the fact that the parallel curved sides of the shell element of Figures 6(a), 6(b) differ by the amount drdϑ, as well as that the different action lines of the forces acting on the remaining two meridian sides are convergent and present a radial component. In the static matrix (67) they appear in the first element of the first row (Ns), in the second element of the first row (N ϑ ), in the third element of the second row (Ts), in the fourth element of the third row (Ms), and in the fifth element of the third row (M ϑ ). It is remarkable that Ns, Ts, Ms are related to the parallel sides, whereas N ϑ , M ϑ are related to the meridian sides.
Whilst the static Equations (67) can be obtained straight-forwardly, the derivation of the kinematic Equations (66), where the terms sinα/r are partially absent (only the terms related to ε ϑ , χ ϑ , and to the meridian actions N ϑ , M ϑ are maintained), is more troublesome. From the Principle of Virtual Work and the static Equations (67), it is possible to derive the kinematic Equations (66). Such demonstration is reported in the Appendix [6].
It is important to emphasize that different terms of the kinematic matrix operator are determined for the first time in the present treatment based on the static-kinematic duality concept [9].

Static vs kinematic adjoint matrices and identity between FEM stiffness matrix and Ritz-Galerkin matrix
We shall define the Finite Element Method on the basis of the Principle of Virtual Work, and we shall show how this is equivalent to the one based on the Principle of Minimum Total Potential Energy. Let the elastic domain V be divided into subdomains Ve, called finite elements of the domain V, and let each element contain m nodal points. Usually in the two-dimensional cases (plane stress or strain conditions, plates or shells, axisymmetrical solids, etc.), the elements are triangular or quadrangular, with the nodes at the vertices, on the sides and, in some cases, inside. In threedimensional cases, the elements are usually tetrahedrons or prisms with quadrangular sides. To each of the nodal points of the element Ve let there correspond a spline, defined on the sole element Ve if the node is internal, also on the adjacent element if the node is on one side, and also on all the other elements to which the node belongs if this coincides with a vertex. To each node k of the element Ve let there then correspond a diagonal matrix made up of the g vectors if g are the degrees of freedom, i.e., the dimension of the displacement vector. These matrices are referred to as shape functions, and have the following properties: Using the Kronecker symbol, we can write more synthetically The displacement vector can be expressed by interpolation, via the shape functions and on the basis of the nodal displacements: In compact form, the displacement vector field defined on the element Ve may be represented as The two vectors on the right-hand side are the vectors of the equivalent nodal forces, and represent the integrated effect of the forces distributed in the domain and on the boundary of the element Ve. Once the local stiffness matrix is calculated, it would be possible to determine the vector of the nodal displacements {δe} on the basis of the local relation (87), only if the forces {p} acting on the boundary of the element were known beforehand and hence the vector of the equivalent forces {pe} were obtained by integration. Whereas, that is, the body forces {F} are a datum of the problem, the forces {p}, which exchange between them the elements at the reciprocal boundaries, are a priori unknown.
To get round this obstacle and, at the same time, to resolve the general problem of the determination of the vector of all the nodal displacements {δ} of the solid, one must add the relation (87), valid for the element Ve, to all the similar relations valid for the other elements of the mesh. In this way, the surface contributions {pe} all cancel out, except for those that do not belong to interfaces between elements, but which belong to the outer boundary. This operation is called assemblage, and involves an expansion of the vectors {δe}, {Fe}, {pe}, from the local dimension (g × m) to the global dimension (g × n), where n is the global number of nodal points of the mesh. The procedure will therefore be to order all the nodes of the mesh of finite elements, so as to be able to insert the nodes of the generic element Ve in the positions that they should have. This may be achieved by premultiplying the vector of the local nodal displacements {δe} by a suitable assemblage matrix [Ae] T , of dimension (g × n) × (g × m), where all the elements are zero, except for (g × m) elements having the value of unity set in the (g × m) different rows to be filled, and corresponding to the (g × m) columns where the integrals extended to the boundaries of the elements cancel out two by two, since the forces that the interfaces of the elements exchange are equal and opposite. It is easy to verify that the vector (97) has equations as its components. Finally, we thus derive the equation which coincides with equation Where The contributions corresponding to the interface between elements cancel each other out. It is easy to verify that the matrix (105) has equations (101) as its elements, and hence the identity between the global stiffness matrix and the Ritz-Galerkin matrix holds: Static-kinematic duality in the shells of revolution: Reasons for the absence of the terms sin α/r from the kinematic matrix operator As clearly shown by Equations (66) and (67), in the shells of revolution, the static and kinematic matrix operators are not exactly one the transpose of the other but for the change of sign in the algebraic terms, as occurs for deformable three-dimensional solids, beams, and cylindri-cal plates. From a mathematical point of view, this peculiarity derives from the application of the Green theorem to a surface element whose area cannot be expressed simply as the product of two differentials. In the present appendix we will prove this fundamental result, i.e. the static-kinematic duality for the shells of revolution loaded symmetrically, where actually the infinitesimal element area is given by rdϑds.
More in detail, we are going to show how the kinematic equations can be derived from the static equations by means of the Principle of Virtual Work. Let us consider a finite portion of a shell of revolution (S being its surface). Without losing generality, we can assume the portion to be bounded by two parallels and by two meridians (together forming its contour C). The Principle of Virtual Work states that the external virtual work equals the internal one, i.e.: (Nsεs + N ϑ ε ϑ + Ts s + Msχs + M ϑ χ ϑ ) r dϑ ds provided that the static system is equilibrated and the kinematic system is congruent. By means of the static Equations (67), the external distributed forces can be substituted in the first integral at the left-hand side, which thus reads: