屈曲文献On Efficient FE-ulations in the Stability Analysis of Thin-walled Structures paper_ww_fg dvi pdf

COMPUTATIONAL MECHANICS WCCM VI in conjunction with APCOM04, Sept. 5-10, 2004, Beijing, China c ?2004 Tsinghua University Press (1,1);(1,1);(1,1) as follows X = 4 ? I=1 NIXID = 4 ? I=1 NIDI. (9) The nodal vectors XI,DIare generated within the mesh input. With (9)2the initial director D is only orthogonal to at the nodes. Accordingly the current confi guration is approximated in the same way x = 4 ? I=1 NIxId = 4 ? I=1 NIdI, (10) where xI= XI+ uIdescribes the current nodal position vector and dIis obtained by an orthogonal transformation dI= RIDI. The rotation tensor RIis a function of the rotation parameters organized in the vector Iand is computed via Rodrigues formula, see 16 in the context of spatial beams. The fi rst and second variation of RI are specifi ed in detail in 16. h A D C B 4 3 2 1 midsurface ( =0)? ? p ? 3 ? ? X 3 X 2 X 1 Figure 1: Quadrilateralhell element The independent fi eld of stress resultants is interpolated as follows h=N N= 1300Nm 00 01300N b 0 001200Ns Nm = Nb = J0 11J 0 11( ) J0 21J 0 21( ) J0 12J 0 12( ) J0 22J 0 22( ) J0 11J 0 12( ) J0 21J 0 22( ) Ns = ? J0 11( ) J0 21( ) J0 12( ) J0 22( ) ? (11) Here, we denote by 12,13second and third order unit matrices, respectively. The vector contains 8 parameters for the constant part and 6 parameters for the varying part of the stress fi eld, respectively. The interpolation of the membrane forces corresponds to the approach in 17, see also 9.The transformation coeffi cients J0 in (11) are the components of the Jacobian matrix J evaluated at the element center, thus J0 = J( = 0, = 0) = X0,( = 0, = 0) t,with= , (12) using an element wise constant local cartesian frame t . The coeffi cients have to be constant in order to fulfi ll the patch test, see e.g. 18,19. The constants and are introduced to obtain decoupled matrices in the below defi ned matrices H,F and denote the coordinates of the center of gravity of the element. The independent strains are approximated in a similar way h= N (13) with 14 parameters in . The interpolation matrix Ncorresponds to Nreplacing the subscripts by , however the coeffi cients of the third row in Nm = Nb contain the factor 2. The mixed element has to fulfi l the patch test for constant stress states, see e.g. 18. As has been shown in 19 the bending patch test when using above interpolation functions for the stresses and strains can be fulfi lled with substitute shear strains according to 14 but not with the shear strains according to (4). The assumed shear strain fi eld according to 14 reads ? 1 2 ? = J1 ? ? where = 1 2(1 ) B + (1 + )D = 1 2(1 ) A + (1 + )C (14) where the covariant shear strains M with M = B,D and L with L = A,C are evaluated at the midside nodes considering (4)3, see Fig. 1. Thus, the virtual shell strains are given by h G = Bv. (15) Here, v denotes the virtual element displacement vector and B the discrete strain displacement operator which can be derived using the variation of the membrane strains and curvatures according to (7) and the variation of the assumed shear strains (14) along with the fi nite element equations (9) and (10). Inserting above interpolations for the displacements, stresses and strains into the linearized stationary condition (8) yields Lg(h,h),h = numel ? e=1 v T e kg0GT 0HF GFT0 v + fi fa fe fs e ,(16) where numel denotes the total number of fi nite shell elements to discretize the problem. Furthermore the following element matrices are introduced kg= ? e kdAfi= ? e BThdA = GT H= ? e NT CNdAfe= ? e NT(W h) dA F= ? e NT NdAfs= ? e NT (h G h) dA G= ? e NT B dA (17) The matrix k is defi ned by the second variation of the shell strains hT G h= vTkv. The vector of the external loads facorresponds to the standard displacement formulation. The matrices H, F and fecan be integrated analytically, since only polynomials of the coordinates and are involved. The other matrices are integrated numerically by a 2 2 Gauss integration. For the linear case an analytical integration of all matrices is possible, see e.g. 19 for a linear plate. Since the stresses and strains are interpolated discontinuously across the element boundaries the parameters and can be eliminated on the element level. Using Lg(h,h),h = 0 and h?= 0 one obtains =FT1(Gv + fs) =F1(H + fe) (18) and thus ke T v=f ke T =GT HG + kg H = F1HFT1 f=GT( + Hfs+ F1fe) fa. (19) The assembly is performed as within the displacement method. The solution of the global system of equations yields the increment of the global displacement vector V. The update of the nodal displacements is performed in a standard way on the system level, V = V + V = + = + , (20) whereas the stress and strain parameters are updated on the element level. For this purpose the matrices which are necessary for the update have to be stored for each element. EXAMPLES The derived element formulation has been implemented in an extended version of the general purpose fi nite element program FEAP, see Zienkiewicz and Taylor 20. 1. Lshaped beam In the fi rst example we discuss the postcritical behaviour of a clamped Lshaped beam, see Fig. 2. The geometrical data are: length l = 240 mm, width b = 30 mm and thickness t = 0.6 mm. An elastic material with E = 71240 N/mm2and = 0.31 has been chosen, see 21. The postcritical path can be reached for the perfect structure using a branch switching procedure or with the introduction of imperfections, for example a small perturbation load P3= P/1000 in thickness direction. P b l b l Fig. 2 Geometry, initial and deformed mesh at displacementu3= 56.1 mm Table 1. Nonlinear buckling loads mesh for one legelementsSimo et.al.9presentEAS-shell12 bottomcentertopcenter 162681.1371.1371.1981.2581.200 3242721.1281.1901.2481.191 64810881.1251.1861.2441.186 converged solution 91.128 0,0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 0102030405060 Displacement u3 mm Load P N Imperfect: P/1000 Perfect Imperfect, large steps Iterat.EAS-shell 12present 19.3256220E+041.0269311E+05 22.8840619E+031.8419593E+03 31.8272735E+039.2095431E+02 49.8608350E+013.1491938E-01 51.4371684E+027.2982781E-02 61.9857848E+013.6491332E-02 75.2637887E+017.7310783E-07 81.4549719E+014.3099687E-09 95.5153927E+01 101.9751270E+00 116.6634091E+01 124.6571811E-02 132.5106273E+00 143.7558851E-05 153.7154558E-06 165.8605419E-09 Fig. 3 Load de.ection diagram and iteration behaviour of residuum Nonlinear buckling loads are presented in Tab. 1 for diff erent FEmeshes. Here the loading position has a signifi cant infl uence on the buckling load. The complete loaddefl ection behaviour is depicted in Fig. 3 using an arclength scheme on the 68 element mesh. The robustness of the present element in comparison to an enhanced assumed strain element (EASshell 12) is demonstrated by the equi- librium iteration for a large step (from u3= 19.540 mm to u3= 35.614 mm), see Fig. 3. Furthermore the initial and the deformed mesh (for a displacement u3= 56.1mm) are plotted in Fig. 2. 2. Hemispherical Shell with a 18hole The hemispherical shell with a 18hole under opposite loads is a standard example in linear and nonlinear shell analysis. A quarter of the shell is modelled with 1616 elements using symmetry conditions, see Fig. 4. The material properties are E = 6.825 107and = 0.3, the radius is R = 10 and the thickness is t = 0.04. The complete load defl ection curve for a 16 16 mesh is presented in Fig. 5. Results for the present element which are nearly identical with the EASshell 12 show a very good agreement with those reported in 9. Starting with F=0 a maximum load step of 40 is possible with the EASshell 12. For this load step the norm of the residual vector within the equilibrium iteration is given in Fig. 5 and again shows the superior behaviour of the new element. It is important to note that the relative large number of 19 iterations occur for a fi nite rotation element along with large rigid body motions and is not a consequence of the enhanced strain formulation. Moreover, the total load of 100 can be calculated using the present element in one load step with 17 iterations. z 2F yx 2F Fig. 4 Hemispherical shell and deformed mesh for F=100 0 10 20 30 40 50 60 70 80 90 100 0,001,002,003,004,005,006,00 Displacement ui Load F u_x Present u_x Simo 9 -u_y Present -u_y Simo 9 Iterat.EAS-shell 12present element 15.6568542E+015.6568542E+01 22.7885600E+062.8374888E+06 34.6613004E+053.3241348E+05 41.9427725E+052.4512080E+04 56.7170299E+042.8536896E+02 62.6142653E+044.5611620E-02 71.3555091E+041.5785771E-08 83.5529025E+03 95.5833012E+03 109.2807935E+02 114.6902795E+03 122.0239489E+02 132.2367207E+03 141.4962903E+01 152.2588811E+03 169.1847138E-01 171.4030970E+01 185.8607442E-04 199.5610236E-06 Fig. 5 Load deflection diagram and comparison of iteration behaviour 3. Steel girder with holes In the last example we discuss the stability behaviour of a beam with holes in a thin web subjected to a vertical load Pzat the center. Considering symmetry half of the structure is depicted in Fig. 6. The web in the range (94x128.5 cm) is modelled as a rigid plate. The following boundary conditions are taken into account: symmetry conditions at x=0, fully clamped at x=128.5 cm. Furthermore no defl ections in y-direction are permitted at x=y=0. E = 21000 kN/cm2and = 0.3 are chosen as elastic material properties. Since the web is very thin the stability behaviour is governed by local buckles around the holes in the web, see Figs. 8 and 9. The load defl ection curves for three selected points, see Fig. 6, calculated with an arclength scheme and a branch switching procedure at the fi rst buckling load, are presented in Fig. 7. Again the developed element allows relatively large load steps in the postcritical range. The defl ections uy of the fi rst buckling mode and in the buckling range are depicted in Figs. 8 and 9. Similar to an experiment the bending of the fl anges can be observed clearly in a side view of the deformed mesh, see Fig. 10. 31 0.2 1 y z x 1 10.5 10 8 8 21.5 25.5 34.5 94 2 10.5 21.5 1 2 3 25.5 cm Fig. 6 Steel girder with holes 0 20 40 60 80 100 120 140 160 180 200 0,000,200,400,600,80 Displacement ui cm Load P kN u_z (Point 1) -u_y (Point 2) u_y (Point 3) Fig. 7 Load deflection curves -5.810E-01 min -4.681E-01 -3.551E-01 -2.422E-01 -1.293E-01 -1.634E-02 9.658E-02 2.095E-01 3.224E-01 4.354E-01 5.483E-01 6.612E-01 7.741E-01 8.871E-01 1.000E+00 max Fig. 8 First buckling mode -6.472E-01 min -5.407E-01 -4.342E-01 -3.277E-01 -2.212E-01 -1.147E-01 -8.186E-03 9.832E-02 2.048E-01 3.113E-01 4.178E-01 5.243E-01 6.309E-01 7.374E-01 8.439E-01 max Fig. 9 Deformed mesh with displacementuyat P=200 kN Fig. 10 Deformed mesh (amplified by factor 10) at P=200 kN and experiment CONCLUSIONS An effi cient quadrilateral shell element for the nonlinear analysis of thin shells is presented. Due to the applied three fi eld variational formulation and appropriate interpolation techniques for the independent mechanical fi elds the convergence behaviour of the mixed hybrid element is superior compared with enhanced strain elements. This has been illustrated by several numerical examples which include bifurcation and postbuckling response. REFERENCES 1 A.E. Green, P. M. Naghdi, On the derivation of shell theories by direct approach, J. Appl. Mech., (1974), 173176. 2E. Reissner, A note on TwoDimensional FiniteDeformation Theories of Shells, Int. J. Non- Linear Mech., 17, (1982), 217221. 3 S. Ahmad, B. M. Irons, O. C. Zienkiewicz, Analysis of thick and thin shell structures by curvedfinite elements, Int. J. Num. Meth. Engng., 2, (1970), 419-451. 4 K. Noor, T. Belytschko, J. C. Simo eds. Analytical and Computational Models of Shells, ASMECEDVol., 3, (1990). 5 G. M. Stanley, K. C. Park, T. J. R. Hughes, Continuum-Based Resultant Shell Elements, in T. J. R. Hughes, E.Hinton eds. Finite Element Methods for Plate and Shell Structures 1: Element Technology, Pineridge Press, SwanseaU.K. (1986). 6 N. Bchter, E. Ramm, Shell theory versus degenerationA comparison in large rotation finiteelement analysis, Int. J. Num. Meth. Engng., 34, (1992), 39-59. 7 K. K. Liu, S. E. Law, D. Lam, T. Belytschko, Resultantstress degenerated shell element, Comp. Meth. Appl. Mech. Engng, 55, (1986), 261-300. 8 J. C. Simo, D. D. Fox, On a stress resultant geometrically exact shell mo del. Part I: Formulationand optimal parametrization, Comp. Meth. Appl. Mech. Engng, 72, (1989), 267-304. 9 J. C. Simo, M. S. Rifai, D. D. Fox, On a stress resultant geometrically exact shell model. Part I I I: Computational asp ects of the nonlinear theory, Comp. Meth. Appl. Mech. Engng, 79, (1990), 2170. 10 W. Wunderlich: Ein verallgemeinertes Variationsverfahren zur vollen oder teilweisen Diskretisierung mehrdimensionaler Elastizittsprobleme, IngenieurCArchiv, 39, (1970), 230-247. 11 P. Betsch, F. Gruttmann, E. Stein, A 4Node Finite Shell Element for the Implementationof General Hyperelastic 3DElasticity at Finite Strains, Comp. Meth. Appl. Mech.