Hybrid stress and analytical functions for analysis of thin plates bending

Hybrid stress and analytical functions for analysis of thin plates bending

Mohammad Rezaiee-Pajand^I; Mohammad Karkon^II

^ICivil Engineering Department, Ferdowsi University of Mashhad, Iran. [e-mail]

^IICivil Engineering Department, Ferdowsi University of Mashhad, Iran. [e-mail]

ABSTRACT

In this paper, two efficient elements for analyzing thin plate bending are proposed. They are a triangular element (THS) and a quadrilateral element (QHS), which have 9 and 12 degrees of freedom, respectively. Formulations of these elements are based on hybrid variational principle and analytical homogeneous solution of thin plate equation. Independent fields in hybrid functional are internal stress and boundary displacement field. The internal stress field has been calculated using analytical homogeneous solution and boundary field is related to the nodal degree of freedoms by the boundary interpolation functions. To calculate these functions, the edges of element are assumed to behave like an Euler-Bernoulli beam. The high accuracy and efficiency of the proposed elements are demonstrated in the severe tests.

Keywords: Finite element, Hybrid stress functional, thin plate, Triangular element, Quadrilateral element.

1 INTRODUCTION

Plate bending structures are widely used in Civil, Mechanic, and Aerospace Engineering. Finite element method is used to solve many complicated engineering problems; among them, one is analyzing plate bending based on Kirchhoff's theory. In this approach, shear deformations of the structure are ignored (Zienkiewicz, 1977). Up to date, many elements have been suggested to solve thin plates bending. Some of them are formulated based on Ritz's method. In this theory, the efficiency of the element is increased according to its degree of freedom, and the related field polynomials, as well. In other words, increasing the degree of freedom creates complicated elements, which are not so good for analysis because of its soaring expenses. During 1965-1975, researchers tried to create high accurate complicated elements. As a result, they created many different elements to solve plates bending (Argyris et al., 1968; Torres et al., 1986; Bath et al., 1989; Martins, 1997; Rezaiee-Pajand and Akhtary, 1998; Dhananjaya et al., 2010). The idea of multi-field formulation, such as hybrid and mixed ones, emerged and developed widely during these years, (Ghali and Chieslar, 1986; Pian, 1995).

Hybrid functional is established by defining one or more master fields on the element interfaces. This approach provides a powerful tool for creating the high-performance elements. In fact, hybrid elements use several independent fields. The main advantage of this kind of formulation is getting rid of the continuity conditions at the boundaries (Felippa, 2009). The first hybrid stress element for plane problems was proposed by Pian in 1964. In this element, the boundary displacement field and internal stress function were selected as master fields (Pian, 1964). Using hybrid stress functional for analyzing plates bending is similar to the plane problems. Many studies have been made on solving plates bending by means of hybrid formulation. By using Hellinger-Reissner functional, Duan and his colleagues created a five node element in order to analyze plates bending with shear effects (Duan, 1999). Furthermore, Miranda and Ubertini, used hybrid stress functional to provide a simple element for analyzing bending plates based on Reissner-Mindlin theory (Miranda and Ubertini, 2006). Spilker and Munir formulated a hybrid stress element for thin bending plates, too (Spilker and Munir, 1980). In another kind of hybrid functional, known as hybrid Trefftz functional, displacement field is used instead of the stress internal field (Choo et al., 2010; Rezaiee-Pajand and Karkon, 2012). In the following lines, the superiorities of finite element method, based on the stress hybrid formulation, compared with the other schemes are given (Felippa, 2009):

1. Relaxed continuity requirements: Due to use of the boundary displacement field in finite element formulation, this condition has been entered in hybrid functional.

2. Better displacement solution: Employing the boundary displacement field increases the accuracy of structural nodal displacement.

In this study, a quadrilateral element (QHS) and triangular element (THS), which have 12 and 9 degrees of freedom, are suggested to analyze thin plates bending. The hybrid functional is written as the total internal energy and interface energy. For internal one, the total potential energy is used. Independent master fields are the internal stress field and the boundary displacement of the element, respectively. The analytical solution of thin plate differential equation is used for internal stress field. On the other hand, the edges of the element are assumed to behave according to Euler–Bernoulli beam theory. Based on these hypotheses, its deflection and torsion fields are calculated. By depicting these fields in the general coordinates, the boundary fields at the edges of element are produced. Using these boundary fields, the powerful triangular and quadrilateral elements are formulated. Accuracy and efficiency of these elements will be proven by implementing numerical analysis of several different problems.

2 HYBRID STRESS FUNCTIONAL

Based on hybrid variational, the total energy of element can be written in two parts: internal and interface energy. In hybrid stress formulation, the internal functional element is written according to its internal stress field. Moreover, its interface energy can be calculated based on its boundary displacement. The total complementary potential energy (TCPE) is used for inside functional. Furthermore, by using of boundary displacement, the interface functional can be introduced. The (TCPE) principle for linear elasto-statics is as follows

On the other hand, the interface energy of element can be written in the below shape

Therefore, hybrid stress-displacement functional has the coming formula (Felippa, 2009)

In these relations, s_ij is the internal stress of element; C_ijkl is the flexibility matrix, which is the inverse of elasticity matrix; n_j is cosine directions' vector; s_ijn_j is boundary tractions; d_i is the boundary displacement field, and W_P represent the potential of applied loads. The latter is calculated by multiplying nodal forces to related displacements. By making stationary of the hybrid functional (3) with respect to the stress and displacement master fields, the stiffness matrix of the proposed elements can be found.

By using the relation (3), finite element formulation for the suggested elements can be achieved. Usually in the plates bending, moment equations are used instead of the real stress. Therefore, equation (3) for the thin plates will change to succeeding form

In this formula, M is the internal moment vector of element, and T is the vector of boundary tractions on the element's edges. E and v are Young's modulus and Poisson's ratio, respectively. Furthermore, t is the thickness of plate. It should be added that the elastic flexibility matrix C for thin plate bending problems has the subsequent formula

In this matrix, E, v and t are Young's modulus, Poisson's ratio and thickness of the plate, respectively.

3 THIN PLATES BENDING FORMULATION

The differential equation of thin plate bending has the succeeding form (Timoshenko and Krieger, 1959)

Here, P(x,y) is the applied load acting on the plate, and D is the stiffness of plate. According to the figure (1), the boundary traction {T} boundary forces and the bending moments around the axis of x and y, as follows

The shears and moments are determined by the coming equations

In order to find the deflection function, the analytical solution of homogeneous differential equation is used. If the equation in its polar coordinates is solved, the later homogenous answer will be obtained (Herrera, 1984)

By assuming r² = x² + y² and z = x +iy, the deflection function, w_j, will be in the following shapes

The first eleven terms of the expression can be found by utilizing k = 0,1,2. As a result, the succeeding solutions are obtained

It is obvious that in homogeneous solution, there is no advantage over x and y. In other words, the element will be rotational invariant. So the internal moments and boundary tractions can be calculated by this way

In these relations, m is the number of deflection functions. The minimum of terms which are selected from homogeneous solution depends on the element's degrees of freedom. Number of necessary terms (and not adequate) to prevent numerical instability and the rank deficiency of the stiffness matrix is acquired through the next condition (Qin, 2005)

In this relationship n and r, are numbers of nodal degrees of freedom of the element under consideration, and discarded rigid body motion terms, respectively. It should be mentioned, by using the modes related to the rigid body motion leads to spurious strain energy in the element, which they should be removed. The number of rigid body motion modes in the plate is equal to three (r = 3). This includes one transitive and two rotational modes in the direction of x and y. Taking these points into consideration; the formula (19) for the plate changes to the following form

4 BOUNDARY DISPLACEMENT FIELD

To achieve boundary interpolation functions, the beam in Figure (2-1) is assumed to be at the element boundary. This beam has two nodes with three degrees of freedom, (w, q_x, q_y) at each node. As it is illustrated in Figure (2-2), the nodal rotational degrees of freedom are projected onto the beam's longitudinal direction and also perpendicular to it. In other words, the beam has two rotational degrees of freedom, q_n1 , and q_n2 , plus two torsional degrees of freedom, q_s1 and q_s2. In addition, there are two deflections, w₁ , and w₂. To utilize these degrees of freedom, a cubically varying displacement field, and also a linear field for torsion are required to be defined.

The next interpolation functions of the nodal displacements of the beam are found in this way

By assuming a linear interpolation for torsion field, the shape functions for the torsional rotation have the coming form

The boundary fields and are calculated by the following relations

On the other hand, the nodal rotations are projected in the x - y plane, as follows

Furthermore, the boundary field's and , will have the subsequent equations

Also according to figure (3) it can be written

Boundary displacement field is determined by multiplying shape functions in the nodal displacements, as follows

Entries of the interpolation matrix are given in the ^{appendix 1}.

5 FINITE ELEMENT FORMULATION

If the formulas (5), (17), (18) and (30) are input in hybrid functional (4), the following equation is produced

The matrix [G] and [F] can be obtained in the next shapes

In this Formula t is the thickness of the element; {D} and {f} are respectively displacement vector and nodal forces. The determination procedure of the nodal equivalent load for QHS and THS elements, are the same as conventional MZC element and BCIZ one, respectively. By making stationary the functional (31) with respect to the stress and displacement degrees of freedom, the succeeding equality will be resulted

According to equation (34), the vector of unknown displacement function will be found as follows

By substituting this equation into relation (35), the stiffness matrix for proposed element can be established in the below form

6 NUMERICAL TESTS

In this section, the efficiency of authors' elements, THS and QHS, is evaluated by using several numerical tests. In all of these benchmark problems, Poison's ratio is assumed to be equal to 0.3. For a comprehensive comparative study, some elements from other researchers are used, which are presented in the following lines

1. Discrete Kirchhoff triangular element (DKT) (Batoz et al., 1980).

2. Discrete Mindlin triangle plate element (DKMT) (Katili, 1993a).

3. The triangular element according to Mindlin's theory (RDKTM) (Chen and Cheung, 2001).

4. The triangular element with C¹ condition and integration by three-point rule (T6/3) (Zienkiewicz and Lefebvre, 1988).

5. 6-node C⁰ triangular plate bending element (AST6) (Sze et al., 1997).

6. Seven nodes triangle element of second order and assumed strain (MITC7) (Bathe et al., 1989).

7. The triangular element according to Reissner-Mindlin's theory (ARS-T9) (Soh et al., 1999).

8. Discrete Kirchhoff quadrilateral element (DKQ) (Batoz and Tahar, 1981).

9. Discrete Mindlin quadrilateral plate element (DKMQ) (Katili, 1993b).

10. 4-node quadrilateral element with 12 DOF (DRMQ) (Sofuoglu and Gedikli, 2007).

11. Unconventional quadrilateral thin plate bending element with 12 DOFs (REC4) (Zienkiewicz and Cheung, 1964).

12. Quadrilateral element with reduced selective integration (Q4-R) (Malkus and Hughes, 1978).

13. The quadrilateral element with reduced integration (Q4BL) (Zienkiewicz et al., 1993).

14. Bathe and Dvorkin's quadrilateral element (MITC4) (Bathe and Dvorkin, 1985).

15. 4-node quadrilateral element based on Reissner–Mindlin theory (ARS-Q12) (Soh et al., 2001).

16. Quadrilateral incompatible element according to Kirchhoff's theory (ACM) (Tocher and Kapur, 1965).

6.1 Effect of the number of analytical functions

The best number of analytical functions for the displacement field is calculated by numerical testing. To find this, a simply supported square plate under uniform distributed load q is considered. The integer number n in THS-n and QHS-n indicates the number of the used functions. Deflection coefficient of the center for the plate with the side a,(w_c/(qa⁴/100D)), is inserted in Table (1). According to these responses, the minimum number of analytical functions, which satisfies the condition of equation (20), yields the highest accuracy results. Considering the outcomes, the triangular element THS with 7 functions and quadrilateral element QHS with 11 functions have the highest level of accuracy for the plate analysis.

6.2 Square plate analysis

In this section, a square clamped and a simply supported plate subjected to a concentrated central load P and also uniform load q, are analyzed. In each case, central deflection coefficient of the plate with side a, is found. The results for simply supported square plate under distributed loading, (w_c/(qa⁴/100D)), are given in Table (2). In addition, Table (3) shows the numerical outcomes for clamped square plate under distributed loading. Furthermore, Table (4) and ⁽⁵⁾ are shown the results, (w_c/(Pa²/D)), for simply supported and clamped square plate under concentrated load, respectively. For each table, percentage of error is exposed in relevant figures. These figures indicate the efficiency and high accuracy of the proposed elements.

Figure 4

Figure 5

Figure 6

Figure 7

It is worth emphasizing that the suggested elements are insensitive to mesh distortion. For evaluating the sensitivity of the presented elements to the mesh distortion, a clamped and simply supported square plate under distributed load of Figure(8) is analyzed in three different meshes, A, B and C. The comparisons of the numerical results with those obtained by other authors are also demonstrated in Tables (6) to (7). According to these tables, the proposed formulations lead to the more accurate response than the other ones.

6.3 Circular plate

In this test, the effects of two distinct loads consist of; a uniform distributed load q, and a concentrated load P acting at the plate center, are investigated. The plate radius is a, and it is solved for both simply and clamped supports with mesh subdivisions shown in Figure (9). Table (8) shows the central deflection coefficient of simply supported circular plate under distributed load (w_c/(qa⁴/10D)). The outcomes of clamped supported structure under distributed load are written in Table (9). Based on these numerical values, the suggested formulations reach the exact answer with a high speed of convergence. Numerical results of the circular plate with simply and clamped supported under concentrated load (w_c/(Pa²/10D)), are inserted in Tables (10) and ⁽¹¹⁾, respectively. The proposed elements yield good accuracy for the circular plate subjected to concentrated load. Error percentages for the elements with different meshes are drawn in the Figures (10), ⁽¹¹⁾, ⁽¹²⁾ and ⁽¹³⁾.

6.4 Razzaque skew plate

Figure (14) shows a 60º skew plate bending, which was originally studied by Razzaque in 1973 (Razzaque, 1973). This structure has simply supported on two opposite edges and free on the other two edges, and it is subjected to a uniformly distributed load. Each side of this plate has length of a. The exact solution was obtained by Razzaque with 16×16 finite difference mesh, as follows

w₀=0.007945qa⁴/D

This structure has been studied with different meshes. Table (12) shows values of the deflection coefficient at the center of Razzaque plates (w_c/(qa⁴/100D)) for the suggested elements. For comparison, other researchers' answers have been recorded in this Table. Figure (15) demonstrates the error percentage of deflection at the center of Razzaque plate for different elements. In this difficult test, a high accuracy of the authors' formulation is quite clear.

6.5 Morley's skew plate

One of the severe problems in the plate analysis is Morley's skew plate. Figure (16) shows this structure. Morley used polar co-ordinates and a least-squares solution procedure to solve an acute skew plate subject to a uniformly distributed load. It is worth emphasizing that this problem poses severe difficulties for numerical methods, since there is a singularity in the bending moments at the obtuse corner. Morley found the following value for the exact central deflection of this problem (Morley, 1963)

w₀ =0.000408 qa⁴/D

This structure will be solved by utilizing the proposed elements. The answers of other researchers are used for better judgment. Table (13) shows the deflection coefficient (w_c/(qa⁴/1000D)) for error percentage of the central displacement is illustrated in Figure (17). This figure reveals that authors' formulations have a high accuracy in Morley's plate analysis.

6.6 Cantilever skew plate

Figure (18) shows the skew cantilever plate. In this structure, the skew cantilever plate subjected to the uniform loading with angles 20, 40 and 60 are analyzed. For the mentioned cause, none of other researchers' elements can converge to the near exact answer. Table (14) to (16) shows the tip deflection coefficient, (w_c/(qa⁴/Et³)) ,for different meshes. These outcomes clearly demonstrate that authors' elements have a high convergence rate in all angles.

7 CONCLUSIONS

In this study, two efficient elements were suggested for solving thin plates bending. The triangular THS and quadrilateral QHS elements have 9 and 12 degrees of freedom, respectively. To formulate these elements, the hybrid stress functional was used. Independent fields in this solution are internal stress and boundary displacement functions, respectively. The internal stress field was obtained by derivation of analytical solutions of the homogeneous part of thin plates governing differential equation. Furthermore, the boundary displacement field is achieved by assuming the element edges behave as a beam. In fact, the beam interpolation functions have been utilized in the presented formulation. On the other hand, the optimal number of analytical functions, which provides the highest accuracy for proposed elements, was calculated from numerical investigations. It is worth emphasizing that using the optimal number of analytical functions, eliminates the zero energy modes and always guarantees the stiffness matrix with full rank. Finally, several numerical tests were done to assess the efficiency of suggested elements. From comprehensive comparison studies of the authors' elements with the famous ones, it was clearly revealed that the presented formulations could offer a significant efficiency and accuracy. It should be added that the accuracy and convergence rate of the element QHS is higher than the element THS.

In the following, the entries of the interpolation matrix (30) for the boundary displacement field are introduced:

APPENDIX II

In order to evaluate the matrix[F], numerical integration technique is used. For this purpose, the Cartesian coordinates, x and y, are written in terms of the natural coordinates x and h. Therefore, the next relationships are held

The interpolation functions, , for the proposed element in natural coordinates, have the coming form:

In this relation, (x_i, h_i) is the nodal coordinate. Thus, the matrix [F] can be calculated as follows

In the last equation, |J|, is the determinant of Jacobian matrix. This matrix has the following formulas

It should be added that the numbers of Gauss's point for calculating the equation (A-4) are 4×4.

APPENDIX III

The matrix {G} is evaluated on the borders of element. For simplicity, the natural coordinate, s, is used. The subsequent equations can be written for individual edge ij

After calculating the matrix {G} for each side, the matrix {G} for the element is obtained.

• Argyris, J.H., Fried I., Scharpf, D.W., (1968). The TUBA family of plate elements for the matrix displacement method. Aeron. J. Roy. Aeron. Soc., 72: 701-709.
• Bathe, K.J., Brezzi, F., Cho, S.W., (1989). The MITC7 and MITC9 plate bending element. Comput. Struct, 32: 797-814.
• Bathe, K.J., Dvorkin, E.H., (1985). A four-node plate bending element based on MindlinReissner plate theory and mixed interpolation. Int. J. Num. Meth. Eng., 21: 367-383.
• Batoz, J.L., Bathe, K.J., Ho, L.W., (1980). A study of three-node triangular plate bending element. Int. J. Num. Meth. Eng., 15: 1771-1812.
• Batoz, J.L., Tahar, M.B., (1982). Evaluation of a new quadrilateral thin plate bending element. Int. J. Num. Meth. Eng., 18: 1655-1677.
• Chen, W., Cheung, Y.K., (2001). Refined 9-Dof triangular Mindlin plate elements. Int. J. Num. Meth. Eng., 51: 1259-1281.
• Choo, Y.S., Choi, N., Lee, B.C., (2010). A new hybrid-Trefftz triangular and quadrilateral plate elements. Applied Mathematical Modeling, 34: 14-23.
• Dhananjaya, H.R., Nagabhushanam, J., Pandey, P.C., Zamin Jumaat Mohd., (2010). New twelve node serendipity quadrilateral plate bending element based on Mindlin-Reissner theory using Integrated Force Method. Structural Engineering and Mechanics, 36(5): 625-642.
• Duan, M., Miyamoto, Y., Iwasaki, S., Deto, H., (1999). 5-node hybrid/mixed "finite element for Reissner-Mindlin plate. Finite Element in Analysis and Design, 33: 167-185.
• Felippa, C.A., (2009). Advanced finite element methods, fourth edition, Mubeen-UET, Lahore.
• Ghali, A., Chieslar, J., (1986). Hybrid finite elements. Journal of Structural Engineering, ASCE, 112: 2478-2493.
• Herrera, I., 1984. Boundary methods: an algebraic theory, Pitman, London.
• Katili, I., (1993a). A new discrete Kirchhoff-Mindlin element based on Mindlin-Reissner plate theory and assumed shear strain fields - Part I: An extended DKT element for thick-plate bending analysis. Int. J. Num. Meth. Eng., 36: 18591883.
• Katili, I., (1993b). A new discrete KirchhoffMindlin element based on Mindlin Reissner plate theory and assumed shear strain fields Part II: An extended DKQ element for thick-plate bending analysis. Int. J. Num. Meth. Eng., 36:18851908.
• Malkus, D.S., Hughes, T.J.R., (1978). Mixed finite element methods-reduced and selective integration techniques: a unification of concepts. Computer Methods in Applied Mechanics and Engineering, 15: 63-81.
• Martins, R.A.F., (1997). A simple and efficient triangular finite element for plate bending. Engineering Computations, 14 (8): 883-900.
• Miranda, S., Ubertini, F., (2006). A simple hybrid stress element for shear deformable plates. Int. J. Num. Meth. Eng., 65: 808-833.
• Morley, L.S.D., (1963). Skew plates and structures, Pergamon Press, Oxford.
• Pian, T.H.H., (1964). Derivation of element stiffness matrices by assumed stress distributions. AIAA J., 2: 1333-1336.
• Pian, T.H.H., (1995). State of the art development of hybrid/mixed finite element method. Finite Element in Analysis and Design, 21: 5-20.
• Qin, Q.H., (2005). Trefftz finite element method and Its applications. Applied Mechanics Reviews, ASME, 58: 316-337.
• Razzaque, A., (1973). Program for triangular bending elements with derivative smoothing. Int. J. Numer. Meth. Eng., 6: 333-345.
• Reissner, E., Stein, M., (1951). Torsion and transverse bending of cantilever plates. TN2369, NACA.
• Rezaiee-Pajand, M., Akhtary, M., (1998). A family of thirteennode plate bending triangular elements. Com. in Num. Meth. Eng., 14: 529-537.
• Rezaiee-Pajand, M., Karkon, M., (2012). Two efficient hybrid-Trefftz elements for plate bending analysis. Latin American Journal of Solids and Structures, 9: 43-67.
• Sofuoglu, H., Gedikli, H., (2007). A refined 5-node plate bending element based on ReissnerMindlin theory. Com. Num. Meth. Eng, 23: 385-403.
• Soh, A.K., Cen, S., Long, Y.Q., Long, Z.F., (2001). A new twelve DOF quadrilateral element for analysis of thick and thin plates. Eur. J. Mech. A/Solids, 20: 299-326.
• Soh, A.K., Long, Z.F., Cen, S. (1999). A new nine DOF triangular element for analysis of thick and thin plates. Computational Mechanics, 24: 408-417.
• Spilker, R.L., Munir, N., (1980). The hybrid-stress model for thin plates. Int. J. Num. Meth. Eng., 15: 1239-1260.
• Sze, K.Y., Zhu, D., Chen, D.P., (1997). Quadratic triangular C0 plate bending element. Int. J. Num. Meth. Eng., 40: 937-951.
• Timoshenko, S., Woinowsky-Krieger S., (1959). Theory of Plates and Shells, McGraw-Hill, New York.
• Tocher, J.L, Kapur, K.K., (1965). Basis of derivation of matrices for direct stiffness method, AIAA Journal 3(6): 1215-1216.
• Torres, J., Amartín, A., Arroyo, V., Díaz del Valle, J., (1986). A C1 finite element family for Kirchhoff plate bending. Int. J. Num. Meth. Eng., 23 (11): 20052029.
• Zienkiewicz, O.C., (1977). The Finite Element Methods in Engineering Science. McGraw-Hill, New York.
• Zienkiewicz, O.C., Cheung, Y.K, (1964). The finite element method for analysis of elastic isotropic and orthotropic slabs. Proc. Inst. Civ. Eng., 28: 471-488.
• Zienkiewicz, O.C., Lefebvre, D., (1988). A robust triangular plate bending element of the ReissnerMindlin type. Int. J. Num. Meth. Eng., 26: 1169-1184.
• Zienkiewicz, O.C., Xu, Z., Ling, F.Z., Samuelsson, A., Wiberg, N.E., (1993). Linked interpolation for Reissner-Mindlin plate element: part Ia simple quadrilateral. Int. J. Num. Meth. Eng., 36: 30433056.

appendix 1

Publication Dates

• Publication in this collection
  30 Jan 2014
• Date of issue
  Aug 2014