Finite element analysis of stiffened plates

The paper presents the development of a new plate/shell stiffener element and the subsequent application in determine frequencies, mode shapes and buckling loads of different stiffened panels. In structural modelling, the plate and the stiffener are treated as separate finite elements where the displacement compatibility transformation takes into account the torsion – flexural coupling in the stiffener and the eccentricity of internal (contact) forces between the beam – plate/shell parts. The model becomes considerably more flexible due to this coupling technique. The development of the stiffener is based on a general beam theory, which includes the constraint torsional warping effect and the second order terms of finite rotations. Numerical tests are presented to demonstrate the importance of torsion warping constraints. As part of the validation of the results, complete shell finite element analyses were made for stiffened plates.


Introduction
Many engineering structures consist of stiffened thin plate and shell elements to improve the strength/weight ratio.The buckling and vibration characteristics of stiffened plates and shells subjected to initial or dead loads are of considerable importance to mechanical and structural engineers.
Among the known solution techniques, the finite element method is certainly the most favourable.Using the technique where stiffeners are modelled by beam finite elements, Jirousek [1] formulated a 4-node isoparametric beam element including transverse shear and Saint-Venant torsion effects.More recent studies on dynamic and buckling problems of stiffened plates and shells are available in [2]- [5].It is a common feature of these finite element based methods that in order to attain displacement continuity, a rigid fictitious link is applied to connect one node in the plate element to the beam node shearing the same section.This approach neglects the out-of-plane warping displacements of the beam section and, in such cases, the usual formulation overestimates the stiffener torsional rigidity.To eliminate this problem Patel et al. in [5] introduced a torsion correction factor.This is analogous to the shear correction factor commonly introduced in the shear deformation beam theory.
The main objects of the present paper is to propose an efficient procedure modelling the connection between the plate/shell and the stiffener, and as part of it the constraint torsion effect in the stiffener.According to Saint-Venant's theory of free torsion, the cross-section does not generally remain plain and the points can move freely in the direction of the beam axis and the angle of torsion changes linearly with a constant rate.If this torsional warping is restricted by external or internal constraints, then the rate of torsion will also change along the beam axis.The theory of constraint torsion was developed by Vlasov [6].Apart from [7], the author could not find any work in the literature involved in the examination of constrained torsion in the stiffener of complex plate/shell structures.However, the effect is obvious, especially in terms of dynamic and stability phenomena when the global characteristics of a structure are investigated, such as frequencies, mode shapes, or critical loads causing a loss of stability.Investigations of stand-alone beam structures proved that an approximate or more accurate modelling of the torsional stiffness, eccentricity, or mass distribution can considerably modify the results.Theoretically -and practically as well, if there is adequate capacity available -beam-type components in complex structures can also be modelled by flat shell, or even spatial finite elements.Consequently, the size of the model and the number of degrees of freedom will change considerably, increasing the time required for calculations and making the interpretation and evaluation of results more difficult.It is a better solution if the properties of components are improved and the ranges of phenomena possible to be modelled are increased at the element level.
As the objective of this paper is to study the effect of constraint torsion and the coupling methods, the shear deformation of the beam is neglected and the formulation of the stiffener is based on the well-known Bernoulli -Vlasov theory.For the finite element analysis, cubic Hermitian polynomials are utilized as the beam shape functions of lateral and torsional displacements.The stiffener element has two nodes with seven degrees of freedom per node.In order to maintain displacement compatibility between the beam and the stiffened element, a special transformation is used, which includes the coupling of torsional and bending rotations and the eccentricity of internal forces between the stiffener and the plate elements.

Beam element
In this work, the basic assumptions are as follows: the beam member is straight and prismatic, the cross-section is not deformed in its own plane but is subjected to torsional warping, rotations are large but strains are small, the material is homogeneous, isotropic and linearly elastic.
Let us have a straight, prismatic beam member with an arbitrary cross-section as it is shown in Fig. 1.The local x axis of the right-handed orthogonal system is parallel to the axis of the beam and passes through the end nodes N 1 and N 2 .The coordinate axes y and z are parallel to the principal axes, marked as r and s.The positions of the centroid C and shear centre S in the plane of each section are given by the relative co-ordinates y N C , z N C and y C S , z C S .The external loads are applied along points P located y S P and z S P from the shear centre.Based on semitangential rotations, the displacement (specifi-cally, the incremental displacement) vector consisting of translational, warping and rotational components is obtained as: where U and U * are the displacements corresponding to the first and second order terms of displacement parameters: (3)

Fig. 2. Notation of displacement parameters and stress resultants
Displacement parameters are defined at the shear centre S as shown in Fig. 2. Accordingly, u, v, and w are the rigid body translations in the x, y and z directions of point S and α, β and γ denote rigid body rotations about the shear centre axes parallel to x, y and z, respectively.The small out-of-plane torsional warping displacement is defined by the ϑ(x) warping parameter and the φ(r,s) warping function normalized with respect to the shear centre.The φ satisfies the following conditions: In these equations the warping function φ and the shear centre location are the same as in the case of free torsion.For thin-walled sections φ = -ω, the sector area co-ordinate (see Vlasov in [6]).When the shear deformation effects are not considered, the Euler-Bernoulli and the Vlasov internal kinematic constraints are adopted as: where the prime denotes differentiation with respect to variable x.The final form of the virtual work principle for the beam structure subjected to initial stresses is expressed as where L , G , Ge are the linear elastic strain energy, the energy change due to initial stress resultants and the potential Per.Pol.Mech.Eng.
energy due to eccentric initial nodal loads, respectively, and W is the work of load increments on incremental displacements.(For further details, apply to [8] - [11].) The first two terms of total potential (6) can be rewritten as: At this point a new displacement parameter, the overall average of U x linear displacement was introduced as: In Eq. ( 7) E and G are the Young's and shear moduli, respectively.The stress resultants in Eq. ( 8) as shown in Fig. 2 are: N the axial force, V r and V s the shear forces acting at the shear centre, M 1 , and M 2 , M 3 are the total twisting moment and bending moments with respect to shear centre, respectively, and B is the bimoment.The stress resultant M W is known as the Wagner effect.With the notation of Fig. 2: Also, sectional properties are defined as It should be noticed that energy functional ( 6) was consistently obtained corresponding to semitangential internal moments because the term (8) due to initial bending and torsion moments was derived based on inclusion of second order terms of semitangential rotations in Eq. ( 3).The detailed derivation of L and G may be referred to -among others - [8], and [10].
The third term of Eq. ( 6) is the incremental work of initial loads: where p and q are the initial surface and volume loads.Considering conservative initial external forces F x , F y and F z acting at material point P (y S P , z S P ) as signed on Fig. 1 of the i-th nodal section, furthermore, assuming that the additional external moments are of semitangential nature, the incremental work of these actions is Definitions of semitangential moments and extensive discussion about their incremental behaviours may be referred to Kim at al. work [8].In [8] the authors justified that the potential energy (8) due to initial stress resultants corresponds to semitangential bending and torsional moments.
For time dependent dynamic problems, volume load increment in the fourth term of Eq. ( 6) is the inertia force and the appropriate external work increment, for beam structure vibrating harmonically with the circular frequency ω, can be written in the following form where ρ is the mass density per unit volume.Substituting the linear displacements from Eq. ( 2) and noting the definition of section properties in Eq. ( 11) and the integral identities of Eq. ( 4), the following expression is obtained for M : (13.b) Finite element analysis of stiffened plates 107 2007 51 2

Finite element model
The derivation of element matrices is based on the assumed displacement field.The nodal vector of seven local displacement parameters is defined as A linear interpolation is adopted for the axial displacement and a cubic Hermitian function for the lateral deflections and the twist: where: Substituting the shape functions into Eqs.( 7), (8), and (13.b) and integrating along the element length L, elementary matrices can be defined as: The explicit -exactly integrated -(14 × 14) element k L linear stiffness and k G geometric stiffness matrices can be found in [10] and the m consistent mass matrix in [12].
Using the lumping technique proposed by Archer at al. in [13] the following lumped, but not diagonal mass matrix can be derived: The numerical accuracy of this rotationally consistent lumped mass matrix was analysed in detail in [12].where

Stiffener transformation
Coupling of the structural components and composition of the system matrices referring to the entire structure are based on the fact that the parameters (degrees of freedom) of connecting nodes are identical.This condition, if it was formulated with the required accuracy, ensures the displacement continuity along connecting surfaces (lines, points).
Majority of publications using the finite element analysis of stiffened panels, where the stiffeners are modelled using beam elements, the beam nodes are forced to undergo the displacements and rotations prescribed by the corresponding plate/shell nodes.In this case the stiffener element follows an edge of the shell/plate element and the constraint condition is introduced by considering a rigid fictious link between the beam section and the plate/shell node N on the common normal.In such a model, on the basis of U displacement vector in Eq. ( 2), the displacements and rotations of a nodal point N, with the co-ordinates r = -y N C and s = -z N C (see Fig. 1) in the plane of the cross-section, will be as follows: where u x , u y , u z , x , y , z are the nodal local displacements and rotations.From the above, the transformation matrix between the local and nodal parameters can be specified as: This transformation takes the eccentricity into account but obviously neglects the effect of torsional warping.

Continuity of rotations
If a beam is connected to another component not only in its cross-section but along a narrow stripe on its surface, the transformation (20) is not sufficient to assure the required displacement continuity.During torsion, while the cross-section turns around point S by an angle α, the originally straight connecting line crossing points N assumes a spiral shape.The rotation arising there is proportional with the distance between points S and N. Using the notations of Fig. 3, the vector of spiral rotation can be described as (21) Supplementing rotations in Eq. ( 19) by this: which yields the modified matrix of transformation between the displacement parameters: where * = y N C z C S − z N C y C S − φ N and φ N is the value of torsional warping function in the node.The only difference between transformations (20) and ( 22) can be seen in column seven.These members link the axial -tensile and bendingmotions with the warping parameter.Accordingly, as regards transformations, a definite difference should be made between beam-to-beam coupling in beam structures when Eq. ( 20) can be applied, and stiffener element coupling, when Eq. ( 22) is suitable.In this form, the latter can be used for any other beam finite elements regardless of the number of element nodes, or the beam theory applied.

Eccentricity of internal forces
The calculation of k Ge load stiffness matrix of the stiffener element requires some remarks.The stiffener load is not known directly as the proportion of the total external (initial) load on the stiffening element depends on relative stiffness conditions.Nevertheless, initial internal forces or contact forces between the stiffener and the plating along the contact line can be calculated from the equilibrium condition of initial state.Hereinafter the contact point should be the node N and using the notation as indicated in Fig. 1, the load eccentricities, if N = P are: ) There is a simple way to calculate the stiffener load stiffness, if the cubic elements are used to define the initial stress state.It follows from the shape functions (15) that the N normal and V r , V s shear forces (see Fig. 3) are constant along a straight beam element, but different from element to element, and the bending and torsional moments are linearly varying.This internal force distribution can be replaced by external forces acting at the two end nodes of an element: With these end loads and (23) eccentricities in Eq. ( 18), the additive stiffness due to off axis contact loads acting along the joint line is expressed as from which the k Ge matrix, likewise Eq. ( 18) with P = N, can be derived in a simple way: Here it is worth to note that the contact forces are part of the internal force system and the internal moments -as it was proved in [8] -are of semitangential nature.For this reason the moment terms are missing in Eq. (25).

Numerical analysis and discussions
With the assembled system matrices the equation of motion of the structure without mechanical load increments can be written as This is the general equation and it can be reduced as a special case to get the governing equation for free vibration and buckling problems as follows: where ω is the natural frequency λ is the critical load parameter.
The goal of the following numerical study is to compare the adequacy of two stiffener coupling transformations detailed in Section 3. First is the usual "rigid lever arm" coupling according to Eq. (20), and the other is the proposed stiffener coupling transformation including the internal force eccentricity in accordance to Eqs. ( 22) and (25).In the following these will be called of BM (beam) and ST (stiffener) coupling, respectively.
These transformations, together with the 7 degrees of freedom per node beam element and a four-noded thick shell element, were implemented in the VEM7 finite element system.The fournoded thick shell element was derived by combining a quadrilateral Mindlin plate element of Bathe and Dvorkin [14] (known as MITC4, mixed interpolation of tensorial components) with a plane-stress element where the contribution of drilling degrees of freedom was taken into account as it was proposed by Allman and Cook [15].

Model definition
A rectangular stiffened panel on Fig. 4 consists of a flat plate with equally spaced longitudinal thin walled T-stiffeners that span between girders.Fig. 5 shows a section of the stiffened plate considered in this investigation.Because of the symmetry in the structure, only a portion of the plate of width b (b = 600 mm) with a T-stiffener centred on the plate strip, was modelled.In the finite element models the rotation about the longitudinal axis (the X axis in Fig. 4) and the lateral displacement were suppressed at all the nodes along the longitudinal edges to simulate the panel continuity (symmetry boundary conditions, x rotation and u y displacement are zero) and the X = 0 and X = L ends of the panel are fixed.The material properties are: E = 2, 010 5 MPa, ν = 0,3 and the mass density: ρ = 8,0 10 −9 N sec 2 /mm 4 .
In order to model the wider range of behaviour of the panel, the plate dimensions and the beam section shape unchanged (b = 600 mm, t = 4 mm), the area of stiffener was scaled in proportion to the web thickness.Using the usual thin walled approximations, the cross sectional properties for the T-beam as the function of t w are: and the non-dimensional plate to beam area ratio parameter is given by In order to verify the BM and ST results of present study a COS-MOS/M shell model was employed.In that model the plate and the thin walled beam was composed of the same four node shell4t thick shell element with six degrees of freedom per node.The mesh of the skin plate was 36x12 for both VEM7 and shell   ).This mesh size found to be sufficient to attain converged results up to three digits.

Dynamic analysis
The mode shapes of the stiffened panel for the first three modes are shown in Fig. 6 and Fig. 7 presents the variation of the first three frequencies in terms of δ size parameter.As it can be seen the frequency of the b1 bending mode is the same in both BM and ST cases for all δ value.In contrast with bending modes, there is a significant difference in the values of BM-t1, ST-t1 and BM-t2, ST-t2 frequencies.The curves on Fig. 7 show that the difference is zero if δ = 0 (no stiffener) and tends to zero with increasing stiffener size and rigidity.In these extreme cases the plate or the stiffener rigidity controls the structure, thus the coupling method is of less importance.It can be stated generally, that ST coupling results in a less rigid model with smaller frequencies.It is observed from the Tables 1.a,b that in case of δ = 0,2 the rate of decrease is around 17%.The COSMOS/M results are marked with black dots on Fig. 7 Finite element analysis of stiffened plates 111 where an excellent agreement has been found for the frequencies between the ST and shell results.

Buckling analysis
To investigate the effect of BM and ST coupling methods on the buckling loads and modes, the elastic buckling of the panel subjected to longitudinal compression is studied in this section.This kind of uniaxial compression can be produced by an U x0 axial displacement of the X= L end of the panel.Here, U x0 = 1 mm compression corresponds to normal stress in the beam and plate, respectively.The buckling mode shapes of the stiffened panel for different stiffener sections are shown in Fig. 8a-b.Fig. 9 shows the change of buckling modes and λ critical load parameter in terms of δ size parameter.If the stiffener section is small the buckling mode is the overall (sometimes called Euler buckling) flexural mode and the result is independent of the coupling method.For higher stiffener sections torsional buckling mode, called of tripping will occur prior to flexural buckling.In contrast with flexural mode, there is a significant difference in the buckling loads of BM-t1 and ST-t1 coupling method, and just as in the dynamic analysis, the ST coupling results in a less rigid model with smaller critical loads.It is observed from the Table 2 that in case of δ = 0,2 the rate of decrease is around 30%.As the stiffener tripping is a coupled lateral torsional-bending motion, the accurate modelling of torsional properties are of great importance.A detailed analysis of the different buckling modes of stiffened panels, including the parametric analysis of tripping, can be found in recent papers of Yuren et al. [16], Sheikh et al. [17] and Hughes et al. [18].With increasing stiffener size and rigidity the difference between BM and ST results vanish.The uniform asymptotic critical load in Fig. 9 indicates the buckling of plate between the stiffeners, as it is shown in Fig. 8. On Fig. 9 quite satisfactory agreement can be seen for the critical load parameter between the ST and shell results marked with black dots.

Conclusions
In this study a detailed numerical evaluation has been performed to prove the efficiency of the proposed stiffener -plate/shell coupling method.It was shown that in all torsion related cases the proposed ST method leads to a less rigid model.The results show good agreement with complete shell solutions.This fact indicates that the application range can be extended.Though further work can be undertaken to perform dynamic and buckling analysis of really curved panels with stiffeners, the newly developed coupling method can be useful for future investigators.

Fig. 7 .
Fig. 7. Change of mode shapes and natural frequencies

Fig. 9 .
Fig. 9. Change of buckling modes and load parameter