Development of the DKMQ Element for Analysis of Composite Laminated Folded Plate Structures

The discrete-Kirchhoff Mindlin quadrilateral (DKMQ) element has recently been developed for analysis of composite laminated plates. This paper presents further development of the DKMQ for analysis of composite laminated folded plates. In this development, a local coordinate system is set up for each element at its centroid. The DKMQ stiffness matrix is superimposed with that of the standard four-node plane stress quadrilateral element to obtain a 24-by-24 folded plate stiffness matrix in the local coordinate system. To avoid singularity of the stiffness matrix, a small stiffness coefficient is added in the entries corresponding to the drilling degrees of freedom. The local stiffness matrix and force vector are then transformed to the global ones and assembled. The accuracy and convergence of the folded plate element are assessed using a number of numerical examples. The results show that the element is accurate and converge well to the reference solutions.


Introduction
Folded plate structures (or, shortly folded plates) are assemblies of flat plates, inclined in different directions, and connected along one or more of their edges.Folded plates are widely used in engineering applications such as stairs, roofing structures, corrugated flooring decks, box girders, ship hulls.Folded plates are often made of isotropic and homogeneous materials such as metals.In line with the advances in material technology, composite laminated materials provide an excellent alternative for designing stiffer, stronger and lighter folded plates.Composite laminated plates are plates formed by stacking layers of different materials [1,2].Published researches on composite materials can be categorized into composite material developments, e.g.[3,4], mechanics, e.g.[1,[5][6][7], and design, e.g.[8,9].This paper belongs to the mechanics" category.
Nowadays the practical method to analyze the mechanical behavior of plates and shells (including folded plates) is the finite element method (FEM).Since the early development of the FEM in the late 1950"s, various plate bending elements have been developed and applied in engineering practice.
The commonly used underlying plate theories in the development of plate bending elements are the Kirchhoff theory and Reissner-Mindlin theory.The latter includes the effects of transverse shear deformation and rotary inertia and thus is a more general approach than the former.
Among countless available plate bending elements, the discrete-Kirchhoff Mindlin quadrilateral (DKMQ) element proposed by Katili [10] is of our interest.The DKMQ has four nodes and three degrees of freedom (d.o.f) at each node, that is, a deflection and two rotations.It is formulated based on the Reissner-Mindlin theory and the assumed shear strain fields and thus applicable for analysis of both thin and thick plates.The benefits of the DKMQ are: (1) it is free from shear locking; (2) it has only three rigid-body modes (no spurious energy modes); and (3) it passes the constant-curvature and constant-shear patch tests.The DKMQ element has been enhanced to a 24-d.o.f shell element that can take into account the warping effects and coupling bending-membrane energy [11].It has also been enhanced for analysis of composite laminated plates [6].Recently, the DKMQ element has been successfully applied to plate bending buckling problems [12].
The applicability of composite laminated materials to folded plate structures has been increased significantly [13][14].The design and analysis of composite laminated folded plates require an enhancement of the computational methods previously directed to plates with isotropic and homogeneous materials.However, literature addressing bending analysis of laminated composite folded plates is very little.For this reason and considering the good performance of the DKMQ element in analyzing composite laminated plates [6], in this paper we extend the application of the DKMQ element to bending analysis of composite laminated folded plates.
The developed DKMQ folded plate element has six degrees of freedom, comprising three translational and three rotational components, at each element node.The stiffness equation is obtained by superposing the equations of the DKMQ element and standard Q4 plane stress element.A small stiffness coefficient corresponding to the drilling degrees is added to the stiffness matrix to avoid equation singularity.

Reissner-Mindlin Folded Plate Theory
Consider a composite laminated folded plate component defined in a global coordinate system X-Y-Z as illustrated in Figure 1.A local coordinate system, x-y-z, is established, where the x-y plane lies on the middle surface of the plate and the z-axis perpendicular to the plate surface.The plate is made from a symmetrically stacking of nl layers of composite materials possessing orthotropic material properties (Figure 2).Each layer k, defined by zk ≤ z ≤ zk+1, is assumed to be in the plane stress condition (the normal stress in z-direction is zero) [2].Applying the Reissner-Mindlin plate theory [1,2,6] and employing the local coordinate system as the reference frame, the displacement of a generic point P(x, y, z) in the plate can be expressed as where u, v, w are the displacement vector components of point P in the x, y, and z directions, respectively, u0, v0, w0 are the displacement vector components of the projected point P0 on the midsurface in the x, y, and z components, respectively, and βx and βy are the rotations of normal line passing through points P0P in the planes parallel to the x-z and y-z planes, respectively.The positive directions of the displacement components and the rotations are shown in Figure 1.An illustration of the displacement and rotation of a normal line in the plane parallel to the x-z plane is given in Figure 2.
The strain components in the folded plate can be split into in-plane strains, εp, due to membrane and bending effects, and transverse shear strains, εs, that is, [2] { } In Eqs. ( 2)-( 4), εx and εy are the extensional strains in the x and y directions, respectively; γxy is the in-plane shear strain; γxz and γyz are the transverse shear strains in the planes parallel to the x-z and y-z planes, respectively.Vectors εm and εb respectively contain the strains corresponding to the membrane and bending deformations.Vector κ contains the curvatures.
Consider a layer k of the composite laminated plate with orthotropic axes L, T, z, and isotropic in the plane T-z (Figure 3).The L-axis is parallel to the longitudinal fiber direction of the composite material whereas the T-axis is transversal to the fiber direction.The relationships between the in-plane stresses and shear stresses in the orthotropic L-T-z axes, that is σL, σT, τLT and τLz, τTz, with their conjugate strains are given as [2,6,15]: where [ ] with (7) In Eqs. ( 7)-( 8), EL and ET are Young"s moduli in the L and T directions, respectively.The Poisson"s ratio ν LT is the ratio of transverse strain in the T direction to the longitudinal strain in the L direction when stressed in the L direction (similarly for ν TL ) [1].The in-plane constitutive matrix, It is apparent that there are totally six independent material parameters, namely, EL, ET, ν LT (or ν TL ), GLT, GLz, and GTz.In many cases, GLz = GLT [2,6] and thus the independent parameters are reduced to five.In a composite laminated plate, the direction of the fiber direction L may change for each layer.This direction is represented by an angle θ between the fiber direction and the local x-axis (Figure 3).For a layer k with the fiber angle θk, the constitutive equations with reference to the local x-y-z axes are [2,6] ; where ; ; Using the definition of the plate stress resultants [1,2,6], Eqs. ( 2)-( 4) and Eq. ( 9), the generalized constitutive equations for the symmetrically laminated plate can be expressed as: where In Eq. ( 14), Nx and Ny are the normal membrane forces in the x and y directions, respectively; Nxy is the shear membrane force; Mx and My are the bending moments corresponding to the normal stresses in the x and y directions, respectively, Mxy is the twisting moment; Qx and Qy are the transverse shear forces on planes normal to the x-and y-axes, respectively.The units of all of these quantities are force or moment per unit width of the plate.In Eqs. ( 15) and ( 16), Epk and Esk are the in-plane and transverse shear constitutive matrices, respectively, as given by Eq. ( 11), for the k-th layer; hk = zk+1 -zk is the thickness of the k-th layer (see Figure 2).In Eq. ( 16), k11, k22, and k12 are shear correction factors to account for the discrepancy between the actual stress distribution and the constant stress predicted by the Reissner-Mindlin theory.These factors depend on lamina properties and lamination scheme [1].
The variational basis of the DKMQ element [6,10] is the Hu-Washizu functional.For the composite laminated folded plate with the middle plate surface A subjected to a static distributed load fz on A (force per unit area), the Hu-Washizu functional can be expressed as: where The independent variables subjected to variations are the generalized displacements, u0, v0, w0, βx, βy, independent (assumed) strains, xz  , xz  , and shear forces, Qx, Qy.

Folded Plate Element Formulation
It is obvious from Eqs. ( 1)-( 3), ( 13), ( 17) and ( 18) that the membrane and bending deformations and stress resultants in the symmetrically laminated folded plate are uncoupled.In view of this, the finite element formulation of the folded plate can be developed simply by superposing a plate bending element with a membrane (plane stress) element.In the present study, the DKMQ element for analysis of composite plate bending structures [6] is superposed with the standard quadrilateral membrane element (Q4) as shown in Figure 4.
The DKMQ element is formulated using the standard bilinear interpolation for displacement component w0 and the bilinear interpolation plus quadratic hierarchical interpolation for the rotations βx and βy [6,10].The assumed shear strain fields are obtained by interpolating assumed tangential shear strains at each element side.These tangential shear strains are obtained using the moment equilibrium and constitutive equations on each side.The shear strain constraint corresponding to the second term of the shear strain energy, Us, in Eq. ( 19) is represented in a discrete manner on the mid-points of the element sides.The detailed formulation leading to the element stiffness matrix of the DKMQ element for a composite laminated plate is presented by Katili et al in [6].
The Q4 element is formulated based on the standard bilinear interpolation of the displacement components u0 and v0.This element corresponds to the membrane strain energy, Eq. (18a), in the Hu-Washizu functional, Eq. ( 17).The formulation leading to the stiffness equations for this standard element can be found in most FEM texts such as [15][16][17].The superposed folded plate element has six degrees of freedom per node in the global coordinate system, that is, U0i, V0i, W0i, θxi, θyi, θzi, i = 1, …, 4.However, it originally has only five degrees of freedom u0i, v0i, w0i, βxi, βyi, per node in the local coordinate system because both the DKMQ and standard Q4 elements have no degree of freedom corresponding to the rotation about local z axis (drilling degree of freedom).To avoid singularity of the global element stiffness matrix, additional degree of freedom βzi per node with a small artificial rotational stiffness has been included (Figure 4).
The procedure to obtain the element stiffness matrix of the present folded plate element is as follows: Step 1: Set up a local coordinate system at the element centroid Position vector of the element"s centroid is (21) in which X1, …, X4 are the position vectors (coordinates) of element nodes 1, 2, 3, and 4, respectively, referred to the global X, Y, Z axes.The local coordinate system is defined by its orthonormal basis vectors î, ĵ, k .Vector î is chosen to be parallel to element edge from node 1 to node 2, that is, where V12 is the vector from node 1 to node 2 and the symbol "| " represents the length of a vector.Vector k , which is perpendicular to the element plane, is where V13 is the vector from node 1 to node 3. Lastly, vector ĵ is Coordinate transformation from a global position vector X to its corresponding local position vector is given as in which Q is the transformation matrix from the global to local coordinates.
Step 2: Calculate the DKMQ and Q4 element stiffness matrices in the local coordinate system The DKMQ element stiffness matrix of the order of 12 × 12 is given as where Bb is the curvature-displacement matrix and Bs is the shear strain-displacement matrix.The full expressions of Bb and Bs are given in Reference [6].
The Q4 element stiffness matrix of the order of 8 × 8 is given as where Bm is the standard membrane straindisplacement matrix (see e.g.[16,17]).
Step 3: Superpose the DKMQ and Q4 element stiffness matrices Each node of a folded plate element has six degrees of freedom (see Fig. 4), that is, three translations u0i, v0i, w0i and three rotations βxi, βyi, βzi, i=1,…, 4. The DKMQ and Q4 element matrices are superposed to obtain the folded plate element stiffness matrix of the order of 24 × 24, kFP.In this process, the entries of the DKMQ and Q4 stiffness matrices are placed on the appropriate location in the folded plate stiffness matrix.The stiffness coefficients corresponding to βzi, i=1, …, 4 are initially zero.
Step 4: Add an artificial stiffness coefficient corresponding to the drilling degrees of freedom Small stiffness coefficients corresponding to βzi are added to kFP.Following Bathe"s suggestion [17], the coefficients are taken to be one-thousandth of the smallest diagonal element of kFP.
Step 5: Transform the element stiffness matrix to the global coordinate system The folded plate element stiffness matrix referred to the global degrees of freedom, KFP, is obtained using the transformation equation where Ti is the transformation matrix from global degrees of freedom U0i, V0i, W0i, θxi, θyi, θzi to local degrees of freedom u0i, v0i, w0i, βxi, βyi, βzi.

Numerical Tests and Discussions
The folded plate element formulation has been coded using Matlab Version 8.3.In this section, a series of patch tests are first presented to verify the correctness of the element formulation and its implementation in Matlab.Then, two folded plate examples taken from Peng et al. [18] are presented to assess the element accuracy and convergence.

Patch Tests using an Inclined Composite Laminated Plate
Patch tests for the developed element were carried out utilizing an inclined composite laminated plate (a typical component of a folded plate) with the mesh shown in Fig. 5.The plate is composed of three layers of orthotropic materials with the lamination scheme of 0°/0°/0° (taken from an example in [6]).
The thicknesses of the skin layers (layer 1 and layer 3) and the core layer (layer 2) are 0.1h and 0.8h, respectively.The material parameters of the skin layers are: EL = 3.4156 MPa, ET = 1.7931MPa, vLT = 0.44, GLT = 1 MPa, GLZ = 0.608 MPa, GTZ = 1.015MPa, k11 = k22 = 0.3521, k12 = 0.The E and G values of the core layer are 10 times lower than those of the skin layers.The tests include the inclined plate under conditions of constant curvatures, constant transverse shear strains, and constant membrane strains.The thickness is h = 100 mm in the constant curvature and constant membrane strain tests.In the constant transverse shear strain test, however, the thickness is taken to be extremely thick h = 1000 m.The reason to take this extreme thickness is that the condition of constant transverse shear strain with zero curvature is only possible when the plate is extremely thick [19].

Constant Curvature Test
The displacements and rotations at the nodes 1, 2, 3, and 4 were prescribed consistent with a constant curvature condition given by (28a) ; (28b, c) The analysis using the DKMQ folded plate element produced the displacements and rotations with maximum relative errors of the order of 10 −10 and 10 −9 , respectively, and the moments with maximum relative errors of the order of 10 −10 .Therefore, the element passes the constant curvature patch test.

Constant Transverse Shear Strain Test
The displacements and rotations at the nodes 1, 2, 3, and 4 were prescribed consistent with a constant transverse shear condition given by ; ; (29) The analysis produced the displacements and rotations with maximum relative errors of the order of 10 −6 and 10 −5 , respectively, and the shear forces with maximum relative errors of the order of 10 −5 .The cause of the errors is the nonzero curvature effect (the curvature will be zero if the plate has infinite thickness) and is not the element inability to produce constant shear stain conditions.Therefore, the element passes the constant shear patch test.

Constant Membrane Strain Test
This test was carried out by consecutively applying a distributed tension force of 0.5 N/m along edge 1-4 in the direction of L-axis, distributed tension force of 0.5 N/m along edge 3-4 in the direction of T-axis, and then a distributed shear force of 0.5 N/m along all edges.In each case, the inclined plate was given an appropriate supports (displacement boundary conditions) to prevent the rigid body motions.The resulting membrane forces in each case were exact.Thus the folded plate element passes the constant membrane strain test.

L-Shaped Folded Plate
An L-shaped folded plate made up of two identical square laminates is subjected to uniformly distributed load fz = −10 Pa on each laminate [18] (see Figure 6 The results using the present element are shown in Table 1.The reference result was scaled from the given figure in Peng et al [18], which was obtained using ANSYS SHELL99 element with 5581 nodes.It is seen that the converged result of the present element is 5.1% lower than the reference solution.A possible reason for this discrepancy is the difference in the plate theory and element formulation used in the present element from those used in the SHELL99 element.Since the discrepancy is relatively small, the present element can thus be regarded to be reasonably accurate.The structure was analyzed using the present element with a mesh of 36 × 36 for each laminate, as shown in Figure 7(b).The deflections of laminate 1 along line X = 0.5 m were observed and presented in Fig. 8, together with those given in Peng et al. [18], which was obtained using ANSYS SHELL99 element with 3025 nodes.It is seen that the deflection curve from the present element is in agreement with the reference curve.

Conclusions
The DKMQ plate element has been developed for analysis of composite laminated folded plate structures.The development was carried out by superposing the DKMQ bending element and the standard Q4 membrane element and adding a small amount of stiffness to the drilling degrees of freedom.
The numerical tests show that the developed element passes all of the patch tests for a folded plate element, gives accurate results, and converges well to the corresponding true solution.Thus, the present folded element can be used as a tool to analyze a composite folded plate in practice.

Figure 1 .
Figure 1.A Composite Laminated Folded Plate Component, the Global and Local Coordinate Systems, and the Sign Convention for the Displacement Components and Rotations

Figure 2 .
Figure 2. Layers (laminas) of Orthotropic Materials in a Composite Laminated Folded Plate and an Illustration of the Displacement and Rotation of a Normal Line

Figure 4 .
Figure 4. Four-node Folded Plate Element Formed by Superposing the DKMQ Element and the Standard Q4 Element (typical degrees of freedom at each node is shown at node 2)

Figure 5 .
Figure 5. Finite Element Model of an Inclined Composite Laminated Plate for Patch Test (a)).The folded plate is clamped along edges a and b.The laminates made up of four layers of orthotropic material with the lamination scheme of −45°/45°/45°/−45°.All of the plies have the same thickness of 0.025 m with the material properties: EL = 2.5 × 10 7 Pa, ET = 1.0 × 10 7 Pa, νLT = 0.25, GLT = GLz = 5 × 10 5 Pa, GTz = 2 × 10 5 Pa.The shear correction factor is taken to be k11 = k22 = k12 = 5/6.The structure was analyzed using different degrees of mesh refinements from 2 × 2 to 64 × 64 meshes for each laminate.The model with 8 × 8 mesh for each laminate is shown in Figure 6(b).The observed result was the deflection at the center point of the top plate.

Figure 6 .
Figure 6.L-shape Folded Plate and its Finite Element Model

1
[18] 4.Half of a Box StructureThis is the second example presented by Peng et al.[18].The problem is a half of a box structure made from three square laminates of the size 1 m × 1 m × 0.05 m as shown in Figure7(a).The box is subjected to a uniformly distributed load fz = −10 Pa on all of its three faces.It is pinned (U=V=W=0) at points A, B, and C. All of the plies in the laminates have the same material properties and shear correction factors as in the first example.We consider the lamination scheme of Case 2 of Peng et al.[18], that is, laminate 1 is taken to be −45°/45°/45°/−45° and laminates 2 and 3 to be 45°/−45°/−45°/45°.(a) A Half Box Structure[18] (b) Finite Element Mesh of

Figure 7 .
Figure 7.A Half Box Folded Plate and its Finite Element Model

Table 1 .
Deflection at the Center Point of the Top Plate (mm) obtained using Different Meshes of the DKMQ Folded Plate Element