Mixed Finite Element Matrix for Plate Bending Analysis Using Weighted-Residual Integral and Weak Formulation

The Conventional Finite Element Method for the analysis of plate bending based on a the principle of stationarity of total potential energy which is well known as Stiffness Formulation , using a variety of conformal and noconformal elements was widely applicable in early of numerical solutions for such problems , later in order to improve the efficiency of such solutions , a more general and flexible formulation called Mixed Formulation based on variational principles which can be regarded as an extension of the Stiffness principle become important alternative. The explicit stiffness matrix for plate bending was given in the most if it`s not in all literatures so the reader can easy follow the solving procedure for such problems and verifying any published results, however such like-matrix (Augmented Matrix) in case of mixed formulation not given, at least for simple elements in the literatures as well as concerned researches. This quietly leads to some difficulties concerning the verification as well as understanding the published results. The main objective of this paper is to introduce this matrix in abbreviated size followed by a applicability through a some detailed examples for some bending models. The derived matrix will be helpful subject of research work in addition to, reveals a very good feeling with understanding and verification the published results in plus to comparing with the analytical solutions of different plate bending problems as the reader can do that.



Abstract-The Conventional Finite Element Method for the analysis of plate bending based on a the principle of stationarity of total potential energy which is well known as Stiffness Formulation , using a variety of conformal and noconformal elements was widely applicable in early of numerical solutions for such problems , later in order to improve the efficiency of such solutions , a more general and flexible formulation called Mixed Formulation based on variational principles which can be regarded as an extension of the Stiffness principle become important alternative. The explicit stiffness matrix for plate bending was given in the most if it`s not in all literatures so the reader can easy follow the solving procedure for such problems and verifying any published results, however such like-matrix (Augmented Matrix) in case of mixed formulation not given, at least for simple elements in the literatures as well as concerned researches. This quietly leads to some difficulties concerning the verification as well as understanding the published results. The main objective of this paper is to introduce this matrix in abbreviated size followed by a applicability through a some detailed examples for some bending models. The derived matrix will be helpful subject of research work in addition to, reveals a very good feeling with understanding and verification the published results in plus to comparing with the analytical solutions of different plate bending problems as the reader can do that.

I. INTRODUCTION
In the Conventional Finite Element Method (CFEM), which is based on the Stiffness Formulation, yields a very accurate displacement for most cases, while the stresses (or moments), which are derived as post-calculation from the displacement field, are not continuous across the boundaries. Although this drawback of (CFEM) may be considered as the major one for solving various engineering problems up to date, and in order to have stress continuity across the element boundary, the moments should be included as independent variables in the formulation. Such an approach is referred to as Mixed Formulation, which is a more general and flexible formulation utilizing variational principles that can be regarded as an extension of the stiffness principle taking significant place in modern field of engineering analysis. Although this Mixed Formulation is a matter of many researches for Finite Element Method as well as other techniques, for plate bending problems, but the explicit forms of the Augmented Matrix were not given for some reasons [1]. This paper describes the construction of this augmented matrix and presents its explicit form through systematic way in order to overcome the size difficulty for finite element plate bending models. In absence of such explicit matrices for all, the correctness and fineness of the published results so obtained remains in doubt due to many interpretations, such as the mathematical expressions of matrix itself, threating of boundary conditions, typographical errors. etc. In this work, some of these doubtfully interpretations were showed and thus contributing in comparing and verifying many of results so obtained in many published research, as well as it will help advanced researchers in this field. Hence it can indicate the reason for this summarized detailed work in order to be simple, but a brief recall revision for many of basic concepts may be necessarily convenient for the reader while some detail of proofs was omitted, for the sake of brevity.
II. REVIEW OF SOME BASIC CONCEPTS OF ENGINEERING AND MATHEMATICAL ANALYSIS In general the finite element method is a technique for constructing approximation functions required in an element-wise application of any variational method, and in order to bring this method in easiness view of understanding, it is necessary to start this review from the method of weighted residuals and weighted-integral statements to arrive at the weak formulation of differential equations, because the weak formulation plays a significance role in the classification of boundary conditions into the so-called natural and essential boundary conditions, which in turn are important in obtaining the interpolation functions as well as the nodal degrees of freedom to be considered in the finite element model. The approximate solution of the Governing Differential Equation (GDEq) such as F( ′ , , ), ∈ [a, b] of a problem consists of assuming and substituting a trial function with undetermined coefficients: , in which will satisfy neither the (GDEq) in its domain nor the boundary conditions on its boundary produce a so-called by CRANDALL [2]; the equation residual or simply an error: R( ) = F( ′ , , ) which was abbreviated by [3] as the Residual R(x) in satisfying such equations, classified by this later reference as domain residual and boundary residual.
The determination of such undetermined coefficients to make residual identical to zero everywhere will yield the exact solution of (GDEq), (i.e. this process sometimes called differential formulation), but due the fact that such process does not always result in the required of linearly independent algebraic equations corresponding to those unknowns. CRANDALL [2] introduced a criteria to insure that there are exactly the same number of equations as there unknowns and make this residual stays close to zero throughout the interval [a,b] which in turn provides a first category of techniques for solve such problems which are: 1-Collocation method developed in [4]: This technique depends on choosing as many locations as there are undetermined coefficients and making the residuals R( ) vanishes at those locations. Using Dirac delta function δ( − ), this method can be expressed in integral form as: ∫ ( − ) ( ) = 0 . 2-Subdomain method [5]: Depending on dividing the desired interval into as many as subdomains as there are undetermined coefficients, by calculating and enforcing the average value of the residual in each subdomain to be zero where mathematically, the average value of any function ( ) in an interval [a,b] defined as (∫ ( ) ) / (∫ ) and this vanishes when the numerator is zero. Hence: ∫ ( ) = 0 . 3-Galerkin`s method [6]: It requires that the weighted averages of the residual over the desired interval should vanish, that is (∫ ( ) ( ) ) / (∫ ( ) ) where ( ) is a weight function and this also vanishes when the numerator is zero, hence ∫ ( ) ( ) = 0 ; which mathematically implies that both of ( ) & ( ) are orthogonal to each other in the interval [a,b]. 4-The method of least squares [7]: This method depends on the minimization technique as it is well known in elementary calculus, the necessary condition for a function ( ) to have a maximum or a minimum (extremum) at a certain point is that the first derivative of the function shall vanish at that point, ( ( )⁄ = 0), while for a function of several variables is that all its first order partial derivatives should vanish thus to minimize the integral of square of the residual over the interval [a,b], casted as: These integral-statements suggest the so-called in general, integral-residual methods in which the weighed-integral statements for certain specified problem provide a new function in integral form of (GDEq) of the problem such as ( ) = ∫ ( , , ′ ) , observing that depends on u(x). These forms are referred to as Functionals which are functions defined by integrals, whose integrands in general are functions of dependent variables and their derivatives that are themselves functions of other parameters such as a position in one, two or three dimensions, time, etc. Thus the functional is a function of functions. For example, of this terminology, the most well-known functional is the total potential energy functional provided in solid mechanics.
In Similar manner as in elementary calculus the extremum of functionals can be obtained through a process called variation of the functional providing the so-called stationary condition of the functional which is the one for which the variation of integral is zero, i.e. I(u) = 0 → ( ) = ∫ ( , , ′ ) = ∫ ( , , ′ ) = 0 , where the symbol δ is called variation operator. The laws for the variation of sums, products, ratios, powers and so forth are completely analogous to differentiation operator. However, it is important to note that the differential of a function df in differential calculus represents a first-order approximation to the change in the function along a particular curve while the variation of the functional δF in calculus of variations is the first-order approximation to the change in the functional from curve to curve. This branch of science is termed variational principles, and for detailed reference on this topic, refer to [8]. In the so-called Rayleigh-Ritz analysis method which is considering as a second group of solving techniques , a functional corresponding to a specified problem utilizes and through a substitution of the trial function ( ) into the functional and invoking the stationarity condition (i.e. variational principle), similar to the method of least squares above. A system of algebraic equations for the parameters c (also called Ritz coefficients) can be obtained, and thus Rayleigh-Ritz method is a general procedure for obtaining approximate solutions of problem expressed in variational form which is different from weighted-residual form.
If the position vector of any point on anticlockwise sense oriented curve in xy-plane is parameterized with arc length s then in Cartesian system is:  Integration by parts with respect to y ; yields: if , 1 = ( ) represent a curve 1 and 2 = ( ) represent C 2 then the first two terms of last the integral can be combined as a single line integral over closed = 1 ∪ 2 oriented in counter-clockwise sense (where = ) and then the first required relation is deduced as:

BENDING ELEMENT
The Mixed Formulation for orthotropic thin plate bending element based on Kirchhoff plate theory can be introduced according to Hellinger-Reissner functional, but here a technique of weighted-residual methods will be followed, which are based on weighted-integral statements of the governing equations which in turn reduced to the so-called Weak Formulations as it is introduced by [3] in the following sub titles.

A. Review of Mixed Formulation Principle
The Weighted-residual method is introduced to provide means for obtaining as many independent relations as there are unknown coefficients in approximate solution when it is substituted in a weighted integral statement of a given governing differential equation. These independent relations can be obtained through two main ways as: The first one by direct substitution of a trial function into a weighted-integral statement of the GDEq. such that the residual in GDEq. be orthogonal to the set of weight functions i.e.∫ ( ) ( ) = 0 , which will provide the required number of equations for determining the coefficients of the trial function. The trial function as well as the weight functions must satisfy the continuity requirements by weighted-integral statement and all types of the specified homogenous boundary conditions of the problem, in which it may require higher -order functions. The second way is exactly as the first one but before substitution is made, the even-order of differentiation of dependent variable in GDEq must be equally distributed among the dependent variable and the weight function. This procedure is well known as a weak form of a weightedintegral statement of GDEq. Hence it is clear that the continuity requirements on trial functions will be reduced, and thus the weak-form name is organized. This weakness of the continuity can be performed using Green-Gauss theorem (generalized integration by parts indicated above), therefore as result of this weakened process the so-called natural boundary conditions will be included in the formulation which in turn the trial functions need not be to satisfying the conditions priory, rather than it should satisfy only the so-called essential condition. The construction of weak formulation can be summarized in the following three steps organized by [3]: 1. Moving all expressions of the governing differential equation to one side and multiplying the entire equation with a weight function and recasting it in a homogenous weighted integral form over its domain. 2. Weaken the continuity by trading the differentiation from dependent variable to the weight function in the governing differential equation using Green-Gauss theorem. 3. Identifying of the primary and secondary variables of the weak form, with the aid of role suggested by [3] that is: after completing weaken process, boundary terms will involve both the weight function and the dependent variable then, the actual dependent variable of the problem which is expressed in the same form as the weight function, is called the Primary variable, while the coefficients of the weight function and its derivatives in the boundary expressions are referred to as the Secondary variables, and their specifications in case of stiffness formulation, represents the Natural Boundary Conditions. and the Essential Boundary Conditions respectively.

B. Review of the Basic Relations and Governing Equations for Bending of Thin Plate
An element × of a plate subjected to uniformly distributed load per unit area ( , ) as shown in the Fig. 2, where the positive sign convention for moments (i.e. stress resultants) and forces per unit length are indicated: And then equations of equilibrium as following [9]: The properties of natural anisotropic as well as structural anisotropy material depending on the direction, and in such material all the strains are coupled to all the stresses [10]. The orthogonally anisotropic is an anisotropic material that has symmetrical three mutually perpendicular planes with respect to its elastic properties [11], such plate material is referred to as the orthotropic plates that has wide practical applications in engineering, where if the principle material axes that are normal to the planes of symmetry of 2Dorthotropy coincide with the x & y coordinate axes then the four elastic constants; young's moduli along the principle material axes , and passion's ratios & that characterize the decrease in a specified-direction during force applied in the other one, since in general the lateral strain is proportional to the longitudinal strain [12], are required for full description of such orthotropic material stress-strain relationships, where due to symmetry the Betti's reciprocal theorem reveals = ; therefore = ⁄ , and for sake of simplifying, the following notations will be used hereinafter: = √ , = √ . The bending rigidities , , twisting rigidity and the shear modulus which characterize changes of angles between principle directions x & y, for plate of thickness are given by [14]: and the relationship between bending , twisting moments and applied load is governed by the following differential equation of equilibrium Substitution of expressions (3) into (4) yields the governing differential equation for the deflection of the orthotropic plate as [14]: where = 1 + 2 , and 1 = =

C. Weak Form Formulation For Bending Of Orthotropic Thin Plate
Solving relations (3) for Finally, this process step yields where from equilibrium relations mentioned earlier that is = + , and application of the same treatment to the rest of equations (6), yields the following set of weakform equations as: where the foregoing remark is now corrected as indicated in equation (10.130) in the Reference [15].
Utilizing an appropriate interpolation function for each of the field variables in the problem, which can be written as:  Remark: ( 12 ) * this term was given with minus sign in eq. (10.136) in Reference [15] , which seems to be also a typographical error. Remark: in eq. (10.136) in Reference [15], the term (− 33 ) * was given with positive sign, which seems to be also a typographical error, and it should be as indicated in part 4 above. [ All these vectors will vanish at inter-element boundaries due to equilibrium requirements at the any two adjacent boundaries of any two successive elements, but at outer edges (i.e. boundaries of the plate), these line integrals must be investigated for different support-type edges. For example, by considering a rectangular plate then:

D. Weak Form Formulation for Bending of Isotropic Thin Plate
The same foregoing detailed explanation can be extended to the isotropic plate which is the one that has the same material property in all directions, those are:

E. The Four Node Serendipity Quadrilateral Elements
Zienkiewicz called the elements whose nodes are only on the their boundaries, as 'Serendip family' elements, by referring to the famous princess of Serendip noted for chance discoveries [16] and the most known one is the 4noded (corner-noded) element, and according to the formulation in this article there exist four degrees of freedom DOF at each node { , , , }. Since for each of the field variables of the problem, this element has one unknown nodal value per node point; totally four unknowns, then mathematically four undetermined constants must be employed in the polynomial expression chosen to represent each filed variable ; for example the field variable is approximated by bi-linear polynomial as: Hence such element is termed as linear element, since the variation of shape functions about a boundary is of that order. In this article derivation of shape functions for the linear 4-noded rectangular element (abbreviated as [R4]) with local numbering system as shown in Fig. 3 is reviewed as:   16×16 , for isotropic can obtained by the same Basic-Matrix[ ] just by using isotropic elastic constants indicated above. As a sample, the shape of Overall Augmented matrix[K G ] which is obtained from the implementation of this procedure where 0 is standing for zero-valued terms while letter v standing for non-zero valued one, is:

IV. IMPLEMENTATIONS OF THE DERIVED BASIC MATRIX AND NUMERICAL RESULTS
The above mixed rectangular element is now used with its Augmented Matrix to solve simple problems of square plates of length L with different types of edge conditions, under uniformly distributed load and concentrated load at the center of the plate and the results were compared with those from a many published literatures as well as with the available exact solutions. Table I shows a comparison of displacements, bending moments at center and twisting moment at corner of all simply-supported plate subjected to uniform loading & concentrated load at the middle, with those parameters calculated by [15].   [15]. Table II shows a comparison of displacements, bending moments at center of all clamped square plate subjected to uniform loading & concentrated load at the middle, with those parameters calculated by [15].  Table III shows a comparison of displacements, bending moment at the center and bending moment at the mid of clamped edge of a square plate with two opposite sides clamped and the other two simply-supported (C-C /x & S-S/y), and subjected to uniform loading & concentrated load at the middle, with those parameters calculated by [15].

Remark:
These results obtained by imposing the boundary conditions: w=0, mx=0 at simply side x=0, and w=0 at clamped side y=0 and mxy=0 at both center lines x=a/2 & y=b/2. But naturally is to impose mxy=0 at clamped side y=0.The given results for 6x6 by Reddy [15] are in fact it should be for 5x5. Table IV contains the values of deflection, bending moment and twisting moment at the specified locations, for a square plate with two opposite sides simply-supported along y-direction and the other two opposite sides are clamped and free sides along x-direction (S-S/y & C-F/x). In this case the results were compared for two chooses of imposing boundary conditions, although the second one, in which mxy not imposing as boundary condition at clamped edge converges in somewhat like the first one but starting with more nearly converged value than the alternative, however utilizing the deflection formula that is given by Timoshenko (set of Eqs. d,e,g,h) [14], with the aid of Maple program package, emphasizes that the twisting moment mxy along the clamped edge must be zero as it`s well known, in other words it should be imposed boundary condition thus the reason for why not imposing in this example as well as in table II above which agreed with Reddy results [15], not understandable for the authors. In this context and for more clarifying this situation, the graphs of mx at clamped edge along x-axis and mxy at free edge are showed in Fig. 4 and Fig. 5, where they are not zero-values there, and the value of the first one found from this work as (-0.035459) for mesh 10x20 is in agreement with the given calculated value (-0.0355) in the table (7.2.6) by Reddy [17]. While mxy at clamped edge are in zero-values.  Table V contains the values of deflection, bending moment and twisting moment at the specified locations, for a square plate with two opposite sides simply supported along y-direction and the other two opposite sides are free sides along x-direction (S-S/y & F-F/x). Finally Table VI contains the values of deflection, bending moment and twisting moment, at the specified locations, for a simply-supported square orthotropic (graphite-epoxy) plate having martial properties as:  Finally, the second author recognizing that the deep search carried out by first author in analytical and numerical formulation which lead to the above results and the following conclusion.

V. CONCLUSION
The introduced step by step derivation of explicit augmented finite element matrix using mixed formulation for plate bending, which in turn leads to a relatively small space and size of such an explicit terms of the matrix when are utilized in either by hand or computer programming calculations.
Applying of such augmented matrix in solving plate problems with a necessary discussion on reasonable treatment of imposing boundary conditions as it should be in accordance with actual problem behavior rather than looking for a more attractive values, reveals very good results when compared with either analytical or numerical solutions of plate bending problems in available literature as the reader can possibly verify. This augmented matrix will provide an easy and direct means and a land mark reference in verifying many obtained results and avoiding its incompatibleness wrong interpretations in the past and future research developments in the area of solving plate bending problems.