Development of Finite Element Procedures for Fluid-Structure Interaction Citation Liu, Wing Kam (1981) Development of Finite Element Procedures for Fluid-Structure Interaction. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/4yzv-q979. https://resolver.caltech.edu/CaltechThesis:08202021-235130206 Abstract In this thesis the development of finite element procedures for fluid-structure interaction problems is presented. The areas upon which attention is focused are: numerical transient algorithms which emphasize implicit-explicit finite element concepts; finite element kinematical descriptions for modelling fluid subdomains in fluidstructure interaction problems; finite element methodology for nearly incompressible fluids and solids, and beam, plate and shell structures based upon theories which include transverse shear deformations; and finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. All these nonlinear methodologies have been integrated into a working finite element computer code. A number of numerical examples are presented to demonstrate the effectiveness of these approaches. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Civil Engineering) Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Civil Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Hughes, Thomas J. R. Group: Earthquake Engineering Research Laboratory Thesis Committee: Unknown, Unknown Defense Date: 27 August 1980 Other Numbering System: Other Numbering System Name Other Numbering System ID EERL Report 80-06 Record Number: CaltechThesis:08202021-235130206 Persistent URL: https://resolver.caltech.edu/CaltechThesis:08202021-235130206 DOI: 10.7907/4yzv-q979 Related URLs: URL URL Type Description https://resolver.caltech.edu/CaltechEERL:1980.EERL-80-06 Related Item Technical Report EERL 80-06 in CaltechAUTHORS Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 14331 Collection: CaltechTHESIS Deposited By: Kathy Johnson Deposited On: 20 Aug 2021 23:57 Last Modified: 20 Aug 2021 23:57 Thesis Files PDF - Final Version See Usage Policy. 7MB Repository Staff Only: item control page

CALIFORNIA INSTITUTE OF TECHNOLOGY EARTHQUAKE ENGINEERING RESEARCH LABORATORY DEVELOPMENT OF FINITE ELEMENT PROCEDURES FOR FLUID-STRUCTURE INTERACTION By Wing Kam Liu EERL 80-06 A Report on Research Conducted Under Grants from the National Science Foundation and the Electric Power Research Institute Pasadena, California August, 1980 This investigation was sponsored by Grants No. PFR79-00712 and PFR77-23687 from the National Science Foundation, Division of Problem-Focused Research Applications and the Electric Power Research Institute under the supervision of Professor Thomas J.R. Hughes. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the author and do not necessarily reflect the views of the National Science Foundation or the Electric Power Research Institute. DEVELOPMENT OF FINITE ELEMENT PROCEDURES FOR FLUID-STRUCTURE INTERACTION Thesis by Wing Kam Liu In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 1981 (Submitted August 27, 1980) To my Parents to to to to to to to to to -iii- ACKNOWLEDGMENTS I wish to extend my deepest gratitude to the following: my advisor: Professor T.J.R. Hughes for his assistance, understand- ing, careful guidance and encouragement in al} aspects of my thesis; my undergraduate advisor: Professor A. B. Schultz for his motivation of my pursuing this degree; my former excellent teacher: Professor T. B. Belytschko for his mo- tivation of my interest in finite element methods; my colleague: Dr. T. K. Zimmermann for his advice and encouragement in the development of part of this thesis; the Caltech SOPS (Society of Professional Students), the staff, and faculty for their help and sharing; the Computing Center of the California Institute of Technology for providing computing services; Ruth Stratton for expending so much care and energy in the typing of the manuscript; the Electric Power Research Institute and the National Science Found- ation for their interest and support of this study; and finally, my Parents whose unending support and encouragement have brought me this far. ABSTRACT In this thesis the development of finite element procedures for fluid-structure interaction problems is presented. The areas upon which attention is focused are: numerical transient algorithms which emphasize implicit-explicit finite element concepts; finite element kinematical descriptions for modelting fluid subdomains in fluid- structure interaction problems; finite element methodology for nearly incompressible fluids and solids, and beam, plate and shell structures based upon theories which include transverse shear deformations; and finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. All these nonlinear methodologies have been integrated into a working finite element com- puter code. A number of numerical examples are presented to demonstrate the effectiveness of these approaches. Chapter TABLE OF CONTENTS 1. INTRODUCTION I. Il. General Remarks and Previous Studies I-A. Introduction I-B. Structural/Solid Mechanics I-C. Fluid Mechanics I-D. Interface Techniques I-E. Transient Analysis Scope of the Present Investigation References for Chapter } 2. FUNDAMENTAL CONTINUUM MECHANICS THEORY AND FINITE ELEMENT FORMULATION II. IV. General Description of Governing Equations I-A. General Theory of Mixed Lagrangian-Eulerian Description I-B. Balance Laws (Mass and Momentum) Formal Statement of the Initial/Boundary-Value Problem II-A. Mixed Description Linearized Equations of Motion Using Lagrangian Description III-A. Constitutive Equation III-B. Linearized Variational Equation Galerkin/Finite Element Formulation and Matrix Equations . IV-A. Generalized Mixed Lagrangian-Eulerian Linearized Variational Equation IV-B. Matrix equations IV-C. Definition of di ke for Incompressible Fluid and Slightly Compressible Fluid Page PF wwnro 10 16 16 7 23 27 27 29 29 34 36 36 37 44 Chapter -vi- References for Chapter 2 3. NUMERICAL SOLUTION TECHNIQUES I. II. III. IV. Numerical Integration of the Variational Equations I-A. Selective Integration I-B. "Upwind" Finite Elements by Modifying the Standard Gauss-Legendre Rules I-C. Mass Matrix Implementation of Lagrangian-Eulerian Finite Elements II-A. Generalized Description Using Isoparametric Finite Elements Numerical Integration of Rate Constitutive Equations III-A. Numerical Algorithms for Rotational Term III-B. Numerical Algorithms for Integration of Equation (IIIA-10) Transient Analysis Techniques for the Variational Equations IV-A. Temporal Integration of the Variational Equations IV-B. Computer Implementation of Equations (IVA-1) through (IVA-29) References for Chapter 3 4, NONLINEAR FINITE ELEMENT ANALYSIS OF SHELLS I. Il. Ill. Introduction Geometric and Kinematic Descriptions II-A. Geometry II-B. Kinematics II-C. Laminar Coordinate Systems II-D. Fiber Coordinate Systems Integration of Constitutive Equations III-A. Numerical Algorithms 50 50 50 58 63 64 65 67 68 70 76 78 81 84 92 92 94 94 99 102 102 107 107 -vii- Chapter Page IV. Element Arrays 109 IV-A. Material Tangent Matrix 109 IV-B. Strain-Displacement Matrix 110 IV-C. Transformation Matrices 111 IV-D. Tangent Stiffness Matrix and Internal 112 Force Vector IV-E. External Force Vector 113 IV-F. Stress Resultants 115 IV-G. Mass Matrix 117 V. Shell Elements 118 V-A. Reduced Integration Lagrangian Elements and 118 Heterosis Elements VI. Sample Problems 121 References for Chapter 4 140 5. NUMERICAL EXAMPLES OF FINITE ELEMENT ANALYSIS OF 144 INCOMPRESSIBLE VISCOUS FLOWS BY THE PENALTY FUNCTION FORMULATION I. Introduction 144 II. Sample Problems 147 TII. Remark on Slightly Compressible Formulation for 192 Viscous Fluid Flow References for Chapter 5 193 6. APPLICATION TO FLUID-STRUCTURE INTERACTION 196 I. Introduction 196 II. Contact/Sliding Element 198 III. Numerical Examples 201 References for Chapter 6 218 7. SUMMARY AND SUGGESTIONS FOR FURTHER DEVELOPMENT 220 References for Chapter 7 223 -]- Chapter ] INTRODUCTION I. General Remarks and Previous Studies I-A. Introduction Significant attention has been devoted to the development of fluid-structure interaction procedures over the last two decades. However, many fluid-structure interaction problems are of such complex- ity that the method of analysis must be numerical in nature. Finite element methods enable the modelling of such complicated problems. The term "fluid-structure interaction” refers to either — (1) fluids contained within structures, or (2) structures surrounded by fluids. Examples of (1) are the seismic response of ground supported, cylindrical liquid storage tanks, aerospace applications such as rocket fuel tanks, and the response of a dam due to sudden acceleration into a contiguous reservoir. Examples of (2) are the transient motion of sub- merged or partially submerged structures, flow-induced vibrations, normal and abnormal operations of light water reactor. and pressurized water reactor cores. The applications of finite element methods to fluid-structure interaction problems have been made by several investigators. One of the related current research areas is seismic analysis and design pro- cedures for ground supported, cylindrical liquid storage tanks [1,2]. Another application is the use of toroidal tanks as pressure-Suppression pools in certain designs of boiling water reactors [3]. Presently used methods are mostly based on linear small motion response and do not take satisfactory account of nonlinear effects. For example, in the response -2- of tanks to strong earthquakes, lift-off, finite-amplitude sloshing, and nonlinear and inelastic response are important. Nonlinear fluid-structure finite element procedures have been proposed by Belytschka et al. [4-9] and Donea et al. [10,1])] among others. However, most of these studies are oriented to particular appli- cations, specifically nuclear reactors, and these procedures are not flexible enough for a wide class of engineering applications. The ability to solve general classes of fluid-structure interac- tion problems involving finite deformations depends upon the ability to solve the corresponding uncoupled fluid and structural problems, and also the ability to "interface" fluid and structural subdomains. Even though the development of general finite element procedures for fluid- structure interactions is obviously too complicated to handle as one thesis topic, the author's goal is to make significant strides in this direction. I-B. Structural/Solid Mechanics The success of finite element methods in structural/solid mechanics is well known. The main thrusts of research herein were directed toward both low- and higher-order, two- and three-dimensional continuum elements with particular reference to nearly incompressible materials [%2-15,19-21] and the development of effective general thick/ thin shell elements based upon theories accounting for transverse shear deformation [16-18]. Attention is also focused on reduced/selective integration and allied concepts which have facilitated the development of simple yet effective "displacement" models for these constrained media applications [12-21]. For general nonlinear material models, the -3- numerical integration of constitutive equations requires special atten- tion since it is a major cost in the finite element calculations. Many new techniques have been proposed [22-29] for this task. For nonlinear shell] analysis this is complicated by the fact that the stress components are referred to a rotating basis,and a zero normal stress constraint needs to be enforced with respect to the normal direction of this basis [18]. I-C. Fluid Mechanics In the fluid mechanics area, finite element research is just coming of age. Although much effort has been exerted in recent years, it is fair to say that transient finite element, Navier-Stokes algorithms have not achieved the speed and versatility of existing finite difference methods [30-32]. The basic problems of fluid mechanics involve non- symmetric “convection” operators. In most physical problems of engineer- ing interest, convection is dominant. Unlike solid and structural mechanics which involve well-understood symmetric, positive, spatial differential operators, nonsymmetric "convection" operators are stil] not fully understood [19.33.34]. Due to the ability to eliminate the incom- pressibility constraints and pressure unknowns [19], and also to the better understanding of the convection operators [33,35-37], it is hoped that the finite element technique, due primarily to its geometric versatility, will soon equal, and eventually surpass the better finite- difference methods. I-D. Interface Techniques The simplest Eulerian-Lagrangian interfacing technique goes under the name of “arbitrary Lagrangian-Eulerian" (ALE) technique in the finite -4- difference literature [38]. Donea et al. [10] and Belytschko and Kennedy [9] generalized the procedure in straightforward fashion to finite ele- ments. Hughes et al. [39] proposed an alternative mixed Lagrangian- Eulerian description in which each degree-of-freedom may be assigned to move at a fraction of the fluid particle velocity. Many problems, such as sloshing and surface waves, may be handled by these procedures. To illustrate the basic idea, consider a cylindrical liquid storage tank, under the action of a strong earthquake. Field observations after strong earthquakes indicate: (1) large amplitude free surface sloshing, which often results in roof damage, (2) nonlinear inelastic tank re~- sponse such as the well-known elephant's foot bulge phenomenon, and (3) lift-off of tank from foundations. In order to model fluid-structure problems of this type, the following description May be employed: Lagrangian for the shell structure, Eulerian for the main core of the fluid, a split or mixed description for the free surface, and an interpolated description between the Lagrangian and Eulerian subdomains (see Fig. I-C.1). In this way the motion of the fluid may be represented without gross distortions of the mesh. I-E. Transient Analysis Basically, there are two general classes of algorithms for tran- sient nonlinear problems: implicit and explicit [40,41]. Implicit algorithms tend to be numerically stable, permitting large time steps, but the cost per step is high and storage requirements tend to increase dramatically with the size of the mesh. On the other hand, explicit algorithms tend to be inexpensive per step, and require less storage than the implicit algorithms, but numerical stability requires that small steps be employed. There are some problems for which implicit pinjy 0 Buruipjuod aunjonays axl, - j;eus JO udijdiuosap ysew udlseing - uDIbuDsbD4 1-O1I ‘614 1 uw Suol|disasap paxity /41|ds — Uudl}duosap UDI 4ajnNg [ —¢ 29D44ns 9844 —o-4 ro uolLdiosap uoibupsHp7 []eys -6- algorithms are very efficient and others for which explicit algorithms are very efficient. However, the direct time integration of the equa- tions resulting from fluid-structure interaction problems presents com- putational difficulties that are not normally encountered in single-field problems (structural and fluid). To overcome these difficulties, methods have been developed in which it is attempted to simultaneously achieve the attributes of both classes of algorithms. Belytschko andMullen [42,43] proposed an implicit-explicit nodal partition in which the explicit nodes are inte- grated first in each time step and the results are used subsequently as "boundary conditions" for the integration of the implicit nodes. Hughes and Liu [44,45] proposed an alternate element-by-element implicit- explicit partition which is amenable to stability and accuracy analysis, and, at the same time, can be simply and concisely implemented. In this approach the mesh is divided into implicit and explicit finite element groups only. The notions of "strong coupling," "weak coupling," “interface elements," explicit and implicit nodes [46,47] are completely avoided. The coupling between the groups is fully accounted for by the standard finite element assembly procedures. The convergence, extension to nonlinear problems, and implicit-explicit finite element techniques for symmetric and nonsymmetric systems have been developed in [19,50,51, 54]. A general partitioned transient analysis procedure is proposed by Park [48], which is amenable to a unified stability analysis and incor- porates the mentioned algorithms aS special cases. Various other improvements in transient algorithms have also been recently achieved [49,52,53,55-57]. The applications of these ideas in different branches of engineering problems has proven to be successful [19, 39,42,53,55]. Il. Scope of the Present Investigation This report deals with the developments of finite element proce- dures in fluid-structure interactions. Both geometrical and material nonlinearities are accounted for. In Chapter 2, the fundamental theories in continuum mechanics are reviewed. A general theory of mixed Lagrangian-Eulerian descriptions will be derived. The definitions of convective velocity and referential or mesh time derivative are given. The balance laws, such as conserva- tion of mass and balances of linear and angular momentum are derived with- in the mixed-Lagrangian-Eulerian concept. The formal statements of the initial/boundary-value problem for the mixed description, Eulerian de- scription and Lagrangian description and their corresponding variational equations will be discussed. A class of rate-type constitutive equations is introduced which represents a wide class of material models. The Galerkin/finite element formulation and matrix equations, as well as the incremental constitituve equations are derived. The elimination of the continuity equation and pressure unknown for a viscous incompressible fluid is achieved within the framework of the penalty function formula- tion. A related idea for slightly compressible fluids is also discussed. In Chapter 3, the numerical solution techniques are discussed. The integrals appearing in the variational equations are carried out by numerical integration. Selective/reduced integration procedures are discussed in detail. "Upwind" finite elements are developed based upon modification of standard Gauss-Legendre quadrature rules and also Petrov-Galerkin techniques. These enable stable finite element -8- approximations to highly convective flows. Consistent, lumped and "row~ summed" mass matrices are discussed. The implementation of Lagrangian- Eulerian finite elements is done by introducing a "Lagrange-Euler parameter." The numerical integration of rate constitutive equations is reviewed, and finally, different transient analysis techniques for the variational equations are presented. In Chapter 4, the nonlinear finite element analysis of shells is discussed. The geometry and kinematic descriptions are first defined. "Laminar" coordinate systems and"fiber'"coordinate systems are then introduced in order to account for the zero-normal stress condition and to eliminate the "torsional"degree of freedom, respectively. The numeri- cal algorithms for integrating the constitutive equations for shells are then discussed. The avoidance of "mesh-locking" in the thin shell limit is achieved by a simple modification of the strain-displacement matrix via selective/reduced integration. Different transformation matrices as well as the tangent stiffness matrix, internal force vector, and mass matrix are defined. Different shell elements, such as the selective/ reduced Lagrangian elements and heterosis elements are also discussed. Numerical examples involving a number of the elements are presented to demonstrate the effectiveness of the proposed shell procedures. Numerical examples of incompressible viscous flows by the penalty function formulation are presented in Chapter 5 to demonstrate the ef- fectiveness and accuracy of the proposed formulation. A free-surface wave generation problem is used to demonstrate the mixed Lagrangian- Eulerian theory in Chapter 6. Also presented is an example of practical interest concerning a nonlinear analysis of an inclined cylindrical liquid storage tank. In addition, the axisymmetric buckling of -9- circular cylinders under different boundary conditions due to a uniformly applied axial load is presented. Finally a dynamic analysis of a liquid storage tank is presented. A summary of the present study and sugges- tions for further developments are discussed in Chapter 7. 10. 11. -1Q- References N. Edwards, "A procedure for the dynamic analysis of thin-walled cylindrical liquid storage tanks," Ph.D. thesis, U. of Michigan, Ann Arbor, Michigan, 1969. M. A. Haroun, "Dynamic analysis of liquid storage tanks," Ph.D. thesis, California Institute of Technology, Pasadena, Ca., 1979. M. Aslan, W. G. Godden, D. J. Scalise, "Sloshing of water in torus pressure-suppression pool of boiling water reactors under earthquake ground motions," Lawrence-Berkeley Laboratory reports #7984 and 6754, NUREG/CR-1082-1083, Oct. 1979. T. Belytschko, J. M. Kennedy, and D. F. Schoeberle, "Quasi-Eulerian finite element for fluid-structure interactions," J. Pressure Vessel Tech. 102, 62-69 (1980). T. Belytschko and U. Schumann, "Fluid-structure interaction in light water reactor systems," preprint. T. Belytschko, ed. '"Fluid-structure interactions in LWR systems," First International Seminar on Fluid-Structure Interaction in LWR Systems (International Congress Center, Berlin (West), Germany, August 20-21, 1979), pp. 29-54. T. Belytschko and R. Mullen, “Two-dimensional fluid structure impact computation," preprint. J. M. Kennedy and T. Belytschko, "Finite element method of above core structures," Argonne National Laboratory, ANL/RAS 79-19, 1979. T. Belytschko and J. M. Kennedy, “Computer models for subassembly simulation," preprint. J. Donea, P. Fasoli~Stella, S. Giuliani, J. P. Halleux, and A. V. Jones, "The computer code EURDYN-1M for transient dynamic fluid-structure interaction,” EUR6751EN. Commission of the European Communities Joint Research Center Ispra Establishment, Italy 1980. J. Donea, P. Fasoli-Stella, and S. Giuliani, “Lagrangian and Eulerian finite element techniques for transient fluid structure interaction problems," Trans. Fourth Internat. Conf. on Structural Mechanics in Reactor Technology, San Francisco, August, 1977. 12. 13. 14, 15. 16. 17. 18. 19. 20. 21. 22. 23. -l1- D. S. Malkus and T.J.R. Hughes, "Mixed finite-element methods-- reduced and selective integration techniques: a unification of concepts," Computer Meth. in Appl. Mech. and Eng. 15, 63-81 (1978). I. Fried, "Finite element analysis of incompressible material by residual energy balancing," Internat. J. Solids and Structures 10, 993-1002 (1974). | J. C. Nagtegaal, D. M. Parks, and J. R. Rice, "On numerically ac- curate finite element solutions in the fully plastic range," Computer Meth. in Appl. Mech. and Eng. 4, 153-178 (1974). D. S. Matkus, "Finite element analysis of incompressible solids,” Ph.D. thesis, Boston University, Boston, Mass., 1975. T.J.R. Hughes, M. Cohen, and M. Haroun, "Reduced and selective inte~ gration techniques in the finite element analysis of plates," Nucl. Eng. and Design 46, 203-222 (1978). T.J.R. Hughes and M. Cohen, "The 'heterosis' finite element for plate bending," Computers and Structures 9, 445-450 (1978). T.J.R. Hughes and W. K. Liu, “Nonlinear finite element analysis of shells: Part I. Three-dimensional shells; Part II. Two-dimensional shelts," Computer Meth. in Appl. Mech. and Eng., in press. T.J.R. Hughes, W. K. Liu, and A. Brooks, "Review of finite element analysis of incompressible viscous flows by the penalty function formulation,” J. Computational Phys. 30, 1-60 (1979). T.J.R. Hughes, “Generalization of selective integration procedures to anisotropic and nonlinear media," Internat. J. Numerical Meth. in Eng., in press. D. S. Malkus, "A finite element displacement model valid for any value of the compressibility," Internat. J. Solids and Structures 12, 731-738 (1976). T.J.R. Hughes and J. H. Prevost, “DIRT II --A nonlinear quasistatic finite element analysis program," California Institute of Technology, 1979, R. D. Krieg, "An efficient numerical method for time independent plasticity," SAND 77-0943, Sandia Labs., Albuquerque, N.M., 1977. 24, 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. -12- J. QO. Hallquist, "NIKE2D: An implicit, finite-deformation, finite- element code for analyzing the static and dynamic response of two- dimensional solids," Lawrence-Livermore Laboratory, Report UCRL-52678, University of California, Livermore, March 3, 1979. R. M. McMeeking and J. R. Rice, "Finite element formulation for problems of large elastic-plastic deformation," Internat. J. Solids and Structures 11, 601-616 (1975). H. L. Schreyer, R. F. Kulak, and J. M. Kramer, "Accurate numerical solutions for elastic and plastic models," J. Pressure Vessel Technol. 101, 226-234 (1979). R. D. Krieg and D. B. Krieg, "Accuracies of numerical solution methods for the elastic perfectly plastic model," ASME J. Pressure Vessel Technol. 99, 510-515 (1977). S. W. Key, "Hondo, a finite element computer program for large de- formation dynamic response of axisymmetric solids," Sandia Labo- ratories, Albuquerque, N.M., SLA~74-0039, April, 1974. T.J.R. Hughes and J. Winget, "Finite rotation effects in numerical integration of rate constutive equations arising in Targe-deforma- tion analysis," Internat. J. Numerical Meth. in Eng., in press, A. A. Amsden and C. W. Hirt, "YAQUI, an arbitrary Lagrangian- Eulerian computer program for fluid flow at all speeds,” Report LA~5100, Los Alamos Scientific Laboratory, Los Alamos, N.M., March, 1973. B. Alder (ed.), "Methods in computational physics," in Advances in Research and Applications 3 (Academic Press, New York, 1964). T. Belytschko (ed.), First International Seminar on Fluid-Structure Interaction in LWR Systems (International Congress Center, Berlin (West), Germany, August 20-21, 1979) pp. 29-54. T.J.R. Hughes (ed.}, Finite Element Methods for Convection Dominated Flows, AMD--vol. 34 (ASME, New York, 1979). T.J.R. Hughes, "Current trends in finite element research, Feature article to appear in Applied Mechanics Reviews. 35. 36. 37. 38. 39. 40. 4). 42. 43. 44, 45. -13- T.J.R. Hughes, "A simple scheme for developing ‘upwind' finite ele- ments," Internat. J. for Numerical Meth. in Eng. 12, 1359-1365 (1978). T.J.R. Hughes and J. Atkinson, "A variational basis for ‘upwind’ finite elements," IUTAM Symposium on Variational Meth. in Mech. of Solids, Northwestern University, Evanston, 111, Sept., 1978. A. Brooks and T.J.R. Hughes, "Streamline-upwind/Petrov-Galerkin methods for advection dominated flows," Third Internat. Conf. on Finite-Element Methods in Fluid Flow, Banff, Canada, 1980. C. W. Hirt, A. A. Amsden, and J. L. Cook, "An arbitrary Lagrangian- Eulerian computing method for all flow speeds," J. Computational Phys. 14, 227-253 (1974). T.J.R. Hughes, W. K. Liu, and T. K. Zimmermann, "Lagrangian- Eulerian finite element formulation for incompressible viscous flows." Presented at the U.S.-Japan Conf. on Interdisciplinary Finite Element Analysis, Cornell University, August 7-11, 1978. 0. C. Zienkiewicz, The Finite-Element Method , 3d edition (McGraw- Hill Book Co., 1977). K-J. Bathe and E. L. Wilson, Numerical Methods in Finite-Element Analysis (Prentice-Hall, Inc., 1976). T. Belytschko and R. Mullen, "Mesh partitions of explicit-implicit time integration," in Formulation and Computational Algorithms in Finite-Element Analysis, K. J. Bathe, J. T. Oden, and W. Wunderlich, eds. (M.I.T. Press, Cambridge, Mass., 1977). T. Belytschko and R. Mullen, “Stability of explicit-implicit mesh partitions in time integration,” internat. J. for Numerical Meth. in Eng. 12, 1575-1586 (1978). T.J.R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient analysis: Stability theory," J. Appl. Mech. 45, 371-374 (1978). T.J.R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient analysis: Implementation and numerical examples," J. 46. 47. 48. 49, 50. a1. 52. 53. 54. -14- Appl. Mech. 45, 375-378 (1978). C. A. Felippa, K. C. Park, and H. C. Yee, “Synthesis of staggered solution procedures for coupled field analysis," 2d International Conf. on Appl. Numerical Modeling, Polytechnical University of Madrid, Madrid, Spain, Sept. 11-15, 1978. K. C. Park, C. A. Felippa, and J. A. DeRuntz, "Stabilization of staggered solution procedures for fluid-structure interaction anal- ysis," in Computational Methods for Fluid-Structure Interaction Problems (T. Belytschko and T. L. Geers, eds.), Applied Mechanics Symposia Series (ASME, New York, 1977). K. C. Park, "Partitioned transient analysis procedures for couplted- field problems," presented at the 2d Conf. on Numerical Meth. in Nonlinear Mech., TICOM, University of Texas, Austin, Tex., March 1979, K. C. Park and P. &. Underwood, "A variable-step central difference method for structural dynamics analysis--Part I: Theoretical aspects," ASME Paper No. 79-PVP-120, presented at the Pressure Vessels and Piping Conf., San Francisco, Ca., June 25-29, 1979. T.J.R. Hughes, K. S. Pister, and R. L. Taylor, "Implicit-explicit finite elements in nonlinear transient analysis," Computer Meth. in Appl. Mech. and Eng. 17/18, 159-182 (1979). T.J.R. Hughes and R. A. Stephenson, "Convergence of implicit- explicit algorithms in nonlinear transient analysis,” International Journal of Engineering Science, in press. H. M. Hilber, T.J.R. Hughes, “Collocation, dissipation, and ‘over~- shoot' for time integration schemes in structural dynamics," Earthquake Eng. Struct. Dyn. 6, 99-118 (1978). T.d.R. Hughes, T. K. Caughey, and W. K. Liu, "Finite element methods for nonlinear elastodynamics which conserve energy," J. Appl. Mech. 45, 366-370 (1978). . T.J.R. Hughes, "Implicit-explicit finite element techniques for sym- metric and nonsymmetric systems," Proc. Conf. on Nonlinear Problems 55. 56. 57. -15- in Mechanics, University of Swansea, Swansea, U.K., Sept. 2-5, 1980. T. Belytschko, H. J. Yen, and R. Mullen, "Mixed methods for time integration," Computer Meth. in Appl. Mech. and Eng. 17/18, 259- 275 (1979). P. G. Underwood and K. C. Park, "A variable-step central difference method for structural dynamics analysis, Part II: Implementation and performance evaluation," ASME Paper No. 79-PVP-121, presented at the Pressure Vessels and Piping Conf., San Francisco, CA., June 25-29, 1979. J. P. Wright, "Mixed time integration schemes," Computers and Structures 10, 235-238 (1979). -16- Chapter 2 FUNDAMENTAL CONTINUUM MECHANICS THEORY AND FINITE ELEMENT FORMULATION I. General Description of Governing Equations "A continuum is an abstraction applied to a large collection of material particles. Continuum mechanics deals with deformations and motions of these particles in space and time, under a variety of ex- ternal conditions such as mechanical and thermal processes. The theory of continuum mechanics serves to establish the mathematical formulation of any physical phenomenon that takes place in a material body. This is tantamount to replacement of the actual body with an idealized ‘mathematical body' in the sense that it is not an atomistic view of matter; rather, it associates material bodies with regions of Euclidean space," K. S. Pister [2]. Equations describing the behavior of a continuum can generally be divided into four major categories: (1) kinematic, (2) kinetic (balanced laws), (3) thermodynamic, and (4) constitutive. Detailed treatments of the subject can be found in many standard texts and ar- ticles [1-6]. However, the two classical descriptions of motion, Lagrangian and Eulerian [3,5,6], are not adequate for many free-surface flows and problems of fluid-solid-structure interaction involving finite deformation. The Lagrangian description is generally not suitable for the fluid undergoing large displacements since the mesh will become highly contorted. In addition, when the Eulerian description is used for the fluid it is not compatible with the large displacement of the structure for which the geometries of the boundary are altered. -17- It is the purpose of this section to develop a more versatile description of motion based on the previous developments by Amsden and Hirt [7], Donea et al. [8], Belytschko and Kennedy [9], Hughes et al. [10], and references therein. I-A. General Theory of Mixed Lagrangian-Eulerian Description n Suppose a body B occupies a region ROR S We describe the "sd motion of B by a mapping ¢: R, x [t,t] +R The image of R, at time t is denoted by Rys and the image of z€eR, is denoted by y, j.e., y = $(z,t) (1A.1) Here RS is thought to be the region occupied at t= ty by the "material particles" which occupy Ry at time t, so we wish to think of Ry as the "spatial" region. In order to describe the mixed formulation, we have to introduce a "reference" region. To this end, let us introduce R sd which is fixed throughout. Its motion is defined by the map- n ping $: R, x [t,>t,] +R 54. we assume that Ry is the image of R, at t Rc under the mapping >, and that y is the image of xe Ry. Therefore, y = O(x,t) (1A.2) By composition we may define a third mapping p by x= ylzst) 3 pest) = o'g(-,t),t) (IA.3) The set-up is illustrated in Figure IA-1. We may obtain the classical kinematic descriptions by specializ- ing the definition of R, and d. The Lagrangian description is obtained = >, whereas if R, = Ry and = 1, we have the if> < by picking R, = R, and ~18- Eulerian description. By appropriate selection of RY and ds more gen- eral kinematic descriptions, which are useful for many free-surface and fluid-solid-structure interaction problems, may be obtained. In order to describe this formulation, itis helpful to introduce coordinate sys- tems. For simplicity, let us restrict ourselves to proper orthogonal coordinate frames and let Xao Yao and Zi5 l<ix< Neg denote the Cartesian components of Xs Vs and Z, respectively. Partial differentia- tion operators will be indicated by the following shorthand notation [10]: ‘js = ; ts 7 ("referential derivatives") 4 x. fixed * 3 3 tt ‘ 7 7 " ae 3 lt eS ("spatial derivatives") oY; ot ys fixed (TA.4) i= = ; . = = ("material derivatives") 24 Z, fixed It follows from (IA.4) that: x, = 0 ; Yep 7 O 5 z. = 0 (1A.5) Keg = Gg 7 265 F 855 (1A.6) where 63 is the Kronecker delta. A summary of useful kinematic definitions, associated with the mappings defined above, is given as follows: udljdiudsep udluejng-uolbunsb07 pexiw | Jo} SbulddowW pud suiowop jo wosboip dIyoWsYdS $|- yI ‘64 NIVWOG con So a qaxi4 (“4y 40 NOILVYNDISNOD IN3YYND) SS” NIVWOG IWILVdS Noun NI ATIVYANSAS ; NivWOG TWIYdS.LVW ~ —_NOILOW NI ATAWY NaS NOILOW ? 4 DILYVd -20- y; = o,(z,t) (motion) (1A.7) Us = ¥i> Za (displacement) (1A.8) o4 u; = Y; = a9, /at (velocity) (1A.9) u; = ¥; = a9, /at? (acceleration) (1A.10) Cr = 39; /22, = 855 446 (deformation gradient) (IA.11) (oy, =O, (xt) (motion) (IA. 12) Us = Ya-X, (displacement) (1A.13) $ { us = yi = ag, /at (velocity) (1A.14) uy = yt = 99 /at® (acceleration) (1A.15) Gar = 96, /0x, = 655 t Uae (deformation gradient) (IA. 16) ‘ Ke = v, (zt) (motion) (1A. 17) We FX Zs (displacement) (1A.18) v ( We = Xe = ap. /at (velocity) (1A.79) Wi. = XK = ay, fate (acceleration) (IA. 20) . . = OW. . = 6..+W. i i A, Gin ay ,/ 825 555 Wig (deformation gradient) (IA.21) Equations (IA.7) - (IA.11) are the classical kinematic relations of the Lagrangian description of a continuum (e.g., Us js the particle dis- placement, Us is the particle velocity, and u is the particle accelera- tion). -?1- Convective velocity and material time derivatives Differentiating (IA.12) with respect to t, holding z fixed, and using (IA.17), yields . _ . Ay from which it follows that a7 %4,553 (TA.23) where -=u.- Ul TA. 24 cj = us - (1A. 24) is called the convective velocity. The material time derivative, which appears in continuum conser- vation laws, is the time derivative with the material coordinates (z) held constant as in equation (I1A.4). We assume, however, if f is a scalar function defined on any one of the regions Rye Ry or Ro» via composition, f can be expressed in terms of the other two functions. Without Toss of generality, if f: Rx[t, ty ]>R, then we can define g: Rx[t,»t,]*R and h: Rx[t,st)]>R (using (IA.2) and (1A.3)), by: g(y.t) ~ g($(x,t),t) = F(x,t) (1A.25) n(z.t) = h(y"(x,t) st) = F(x.) In order to avoid introducing new variables (i.e., instead of writing (1A.25)) we will simply write ey = Fly.t) = F(z, t) (1A. 26) referential spatial material or or Eulerian Lagrangian -22- with the understanding that the second and third f's in (IA.26) are ac- tually g and h, respectively. Differentiating (IA.26) with respect to t and holding z fixed, we get fl+fiw, = fy +f, = f ii sti i (IA.27) referential spatial material or or Eulerian Lagrangian Since Ry is fixed throughout in this formulation, it is convenient to express the material time derivative (f) in referential form. Combining (IA.27) and (IA.23), and by the chain rule, we have the important relation fo = fl4¢f 5c (TA. 28) 1 Strain measures Let a = (l-a)z+ay, where 0<a< 1. Also, let ds, ds, be the differential distances between points in RS and Ry respectively. We can write ds? = dz+dz (1A.29) ds, = dy + dy (IA. 30) Consider: z=2(a) and y = y(a) and define: dso - ds’ = 2dal+ E+ da (1A.31) -23- ~ ~ ~ Letting dz = Fda and dy = Fyda. where i = 9z/da, Vy = oy/ 0a» we have: . ) + (4+ (1-20) GE GDI (14.32) ~ ~ {( ni E = | a She Ifa=0, E is called Green or Lagrangian strain tensor. ~ 1, E is called Aimansi or Eulerian strain tensor. ~ If a If the motion of B is rigid, then this is equivalent to saying —E=0 [5]. Also, for a=1/2 it was shown by Hallquist [14] and by Hughes and Winget [13] that E has certain useful numerical calculation proper- ties as will be discussed later. I-B. Balance Laws (Mass and Momentum) Mass density, conservation of mass (equation of continuity) We assume there exists a scalar-valued function p(+,t) called mass density defined on Ry and t)<t<t, with the properties [5]: (1) lyst) > 0 V ye R, (toststy) (ii) If Pc R, isa bounded subregion of R,, and if Py = (Pst) » then m(p,) = | plyst) ay Py is the mass of material occupying Py at time t. We call Py (x) the initial mass density, viz. Pg(x) = p($(x,t,) st) -24- (iii) If P_¢R, is a bounded subregion of R,, and if PY = o(p,.t), then m(P)) = 5 pfy,t) dP, = a p(y.t) Jy dP, y Zz where Jvz = det(dy/az). The principle of conservation of mass states that mass cannot be created nor destroyed if we follow the motion of the given mass [12], j.e., d da _ a m(Py) = $5 | lyst) dPy = 0 (1B.1) y Equivalently: [oJy,]= 0 (1B.2) Carrying out the material time derivative, (IB.2) becomes ) Jz +9 Jy = 0 (IB.3) As Jy2 = Jy, Ue k [5], conservation of mass says: p + p(t) ,-= 0 (1B.4) If the material is incompressible so that the density of each material particle remains constant as it moves, the continuity equation takes the simpler form UK = 0 (1B.5) This is the condition of incompressibility, which has been important in classical hydrodynamics and in plasticity theories. -25- Body forces and contact forces, surface tractions [5] Body forces We assume there exists a vector-valued function b(+,t) (body- force density per unit current volume) defined on Ry such that the resultant force and resultant moment about the origin, 0, of an eqn dimensional Euclidean space due to the body force acting on the mate- rial in Ry are given by | b(y,t) dP 3 | yxb(y,t) dP pew y pl ~~ y y y Here x is a "cross product" or “vector product.” Contact forces We assume there exists a vector-valued function T(+,+,t) (surface- force density per unit current area, actual surface traction, or Cauchy- traction) defined on Ry x U, where U is the set of all unit vectors, such that the resultant force and resultant moment about 0 due to the surface forces acting on the boundary of Pye denoted ape are given by: | Tlysn(y,t),t) dA, s | y X Tlysn(y,t).t) day P p ory ory where n(y,t) is the unit outer normal vector at y € OP. The principles of linear and angular momentum (postulates) [12] Newton's law states that the time rate of change of linear momentum of a given mass is equal to the force exerted upon the mass, i.e., | bly.t)dP,, + | Tly n.t)dA, = a | plyst) aly st)dP, (IB.6) Py apy Py -26- Also, conservation of moment of linear momentum states that the time rate of change of angular momentum of a given mass with respect to a given point is equal to the applied torque referred to the same point, 1.¢e., _ yon, t)dA = ae | yxplyst) ulyst) dP, (1B.7) ~ [ yxbly.tary * Jy xt P Py oP, y Before deriving the mixed Lagrangian-Eulerian equation of motion we have to introduce a very important lemma. Lemma _ IB-1 Let o: R,x [t,t] +R Sd be sufficiently regular. Let Py and p be defined as previously. Suppose g is vector-valued, and g{yst) = g(o(z,t),t) = 9(z,t) then Q|> cr | ply.t) gly,t) dPy = y y [ lyst) a(y.t) a y P Following [5], we can define a second order tensor t{y,t) called Cauchy, or actual, stress-tensor field on Ry: tT, (y>t) are the Carte- sian components of the two-tensor t(y,t). If n(y,t) is the unit outer normal vector of apy at y, then the surface traction vector T is related to the Cauchy stress-tensor T by: T; = TH ni (1B.8) Substituting (IB.8) into (1B.6), using lemma IB-1, and invoking the divergence theorem, we get -27- .. «tb. - pl. = R | (tia b, - pli) dP 0 VPC RY y Continuity of integrand implies 154,57 Dy = oy on Ry ¥ te[t,.t,] (7B.9) Using (IA.28) on the right-hand side term of (1B.9) yields the equations of motion in the mixed Lagrangian-Eulerian description. -s oe .= 7 + iP -C. T5a,gT Pas ol(uy) #4) 505 on Ry i] t<[t,.t,] (1C.10) Using (IB.7) and following [5], we have the symmetry of the Cauchy stress tensor, 1.€., "53 = Ti" Il. Formal Statement of the Initial/Boundary-Value Problem II-A. Mixed Description We will restrict ourselves to fluid flow problems if the mixed description is being used. The object of the initial/boundary-value problem is to find func- tions u,: R, x [t,»t,] +R, u,: RY x[tysty] > R, and p: Rx [t,t +R (pressure) (where R is the closure of R), such that pul, = T4454 + b. (II1A.1) (us) = 0 on Rx [t, >ty] (IIA. 2) Tij . -pd,, + ulus) 5+ (U5) 4] (IIA.3) - g Us = 93 on aRy x [t»t,] (IIA.4) -28- = h i 744M h. on aR x [ty ty] (IIA.5) Us = UG; on Rix [tJ (IIA.6) te = a given representation on Rx [tty] (I1A.7) depending on $, as will be discussed in the next chapter n~ Us = Uo; on Ry (IIA.8) where 9 > 0 is the density, b. is the prescribed body force (per unit volume Ry) is the dynamic viscosity; oF and h. are prescribed bound- ary velocities and tractions, respectively; n. is the unit outward normal vector to aR 3 and Uo and Qo; are the given initial data. i In order to eliminate the pressure unknown and the continuity equation, we employ a penalty-function formulation of the incompressi- bility condition in which (IIA.2) is dropped and the pressure is deter- mined from p = -A(us) (IIA.9) where } > 0 is the penalty parameter; see [15] for a discussion. An alternative is to use a"slightly compressible formulation” in which it is assumed that 0 = p tun + UL aor + puK (IIA. 10) Then by using the equation of state for an isothermal process, F(p,o) = 0 it follows that = 9p _ po Pit 30 ?yt 30 ek yk -29- We further assume 3p/do ~ constant = 8. Then Pit ~ “Bou k and the constitutive equation (IIA.3) becomes t - ’ - ue) .t(us) ; A. [ soi, dt) 65+ wld) 5#CG),¢2 (TAT) 0 where Po is the pressure at t=0. TR “(Py With the continuity equation eliminated, the variational equation corresponding to (IIA.1) through (I1A.8) is | (pliju; +7; 504; dR, = | bju,dR, + I, hju;dA, (11A.12) Ry Ry ary in which — - g F 0 on aRy (IIA.13) By virtue of RY being fixed throughout, we recall U, = (us) + (ui) 56; III. Linearized Equations of Motion Using Lagrangian Description III-A. Constitutive Equation In many engineering applications. it is found convenient to put constitutive equations into rate form. A wide class of plastic, visco- elastic, and viscoplastic models can be put into this form. Nonlinear elasticity can be put into this form by time- differentiating the more usual forms. We will restrict our attention here to linear elasticity, the Saint-Venant-Kirchhoff nonlinear elastic ~30- model, and the Krieg-Key combined isotropic/kinematic hardening plas- ticity model [21]. A class of rate-type constitutive equations for inviscid materials For discussion purposes, consider rate-type constitutive equations of the following form: Ty 7 74g 7 Sague4tesey ~ FyKe(k,2) (IITA.1) where ij = Ci skoY(k 2) (IIIA. 2) Ci ike = Sp jug htoF) (ITIA.3) Signe = (TaeSsk * TyeSik > THKdge 7 THRSig) / 2 (ITIA.4) UK, 2) = (uy + Uy )/ 2 (IIIA.5) Uk, 2] = (Uy g - Up 4/2 (IIIA.6) Here Ci akg is the material constitutive tensor which is a function of Cauchy stress t and deformation gradient F. Uc, 2) is the symmetric part of the velocity gradient, whereas Ur 4] is the skew-symmetric part. Without loss of generality, we further assume the Ciike § possess the minor symmetries: Ci 5ke = “jake = “igek (IITA.7) The tensors Si 5k2 and Ci ike possess the following symmetries: Sigjke ~ Sjike = “Sijek (ITIA.8) * _ ak _ ok Crake = Fike = “Fjek (IIIA.9) The second term on the right-hand side of equation (IIIA.1) {i.e., the -3}- ¥ : it : . ixed 54 5jkeU[k,2] term) is the “rotational part" of the rate and is fixe once and for all by the so-called “objectivity” or "material frame indifference" as it is now often called [22,23]. As has been discussed in the paper by Hughes and Winget [13], there are several ways to account for this rotational term numerically; for details see [13] and [16]. The third term on the right-hand side of equation (IIIA.1) (i.e., the CF sk 2) term) is required to be an "objective tensor," but is otherwise arbitrary. We find it is convenient to employ the Truesdell rate of T7i4 in which Ci akg has the following form: x def Crake = 7 THSket (THe Syk* THeStk * THKSG et THR OG) 2 — (ETTA.10) For further discussion of the "objective tensor," CF ike? the reader is urged to consult [17,24]. Combining (IIIA.1), (IIIA.2) and (IIIA.10), we can write the con- stitutive equation as follows: TH5> Sgke(kee) * SijkeYlkso] (II1A.11) in which ny _ * Crake = Sigkn + FjKe (IIIA.12) Examples of C5 ike, for particular materials under consideration 1. Linear elasticity A wide class of engineering problems can be studied by linear elastic material behavior. However, if finite rotations are involved, the analyst must generalize the smali strain linear elastic model to account for these effects. In this way (IIIA.11) will yield a symmetric -32- stiffness matrix due to the choice of Truesdell rate. The choice CF eg = 0 (Jaumann rate) results in a nonsymmetric stiffness. As has been mentioned in the paper by McMeeking and Rice [19], the CF ike’ will have a non-negligible effect on material response in the high- stress, finite deformation regime. Use of this model in this regime is clearly speculative. If we consider the isotropic linear elastic model, the material constitutive tensor takes the following form: 5. ,82,) (IIIA.13) Crap. = AS;-6 6 u°5k ijkz ij ke + u(s ik° jn * where X,u are the Lamé parameters. 2. Saint Venant-Kirchhoff model The second (symmetric) Piola-Kirchhoff stress tensor Pag is re- lated to the Lagrangian strain tensor Ens by the following relation: [rs 6 tu(6.6 +6 6 JIE Pog pq rs pr-qs ps qr rs = Doars Ens (IIIA.14) where A. are the Lamé parameters. The Ciake are defined by [18,20]: Teo ofoF Fp (IIIA.15) Ci ike =J ip jq kr &s pars where Fay are the Cartesian components of the deformation gradient ten- sor and J is the determinant of F. Carrying out the algebra and de- fining b = FFT (the Finger tensor), (IIIA.15) becomes: _yl Ca agg = FLA by y + u(D; bj 9 + Bs pds) (IITA. 16) ~ 33- 3. Krieg-Key combined isotropic-kinematic hardening plasticity model [21] Using the assumption made in the paper by Krieg and Key [21], we . ' . can write the Ci ike s as follows: = dO. Ci ska, 1g kat MCSik Syg* § 6.,) - 2uo Ns Myo (IIIA.17) jk i2 where A, are Lame parameters. c is set to one during plastic loading, that is, when ¢ > 0; otherwise, < is set to zero. 9 is the Von Mises yield function defined by: = F ((Eq47Ep9)° + (Eqg-Eqg)° + (Eg978 4) + Bot Ehyt Bag ké (IIIA.18) in which Eo = Ta (IITA.19) rs rs rs and the tensor Gag is the location of the center of the Von Mises yield surface, and k is its radius. The plastic hardening rules are: x _2 “Pp and kx = | epg? (IIIA.21) v3 with the effective plastic stretching, dP defined by P _ 2 -P ~P 1/2 d | ie s)Hr,s)! (IIIA. 22) An associated flow rule is used to give: Wess) = EN (ITIA.23) ~34- where Eos = Ens - Bg Ekk / 3 is the deviatoric part of Ens? u z! - (rss) rs (IIIA. 24) 2k“(1+B/3n) and Nas = Els / f2k2(1 +B/3u)7!/2 (IIIA. 25) 8 is a kinematic variable. If 8 is equal to one, we have the case of isotropic hardening; if 8 is equal to zero, kinematic hardening , and if B is in between, it is combined isotropic-kinematic hardening. B is the plastic modulus and it is related to the elastic modulus, E, and plastic modulus, Ess in the uniaxial true stress-strain curve by: _ 6B | Ey = Ey (IIIA. 26) III-B. Linearized Variational Equation Before we can derive the appropriate form of the variational equation, we have to write the constitutive equation in incremental form. Express t as At/At and u as Au/At, where At is the time incre- ment. Equation (IIIA.11) becomes: Following techniques in [17-20] for deriving linearized varia- tional equations, the linearized version of equation (IIA.27) is given as follows: -35- | pAU.U. dR, \ R y R y = | b.u. dR, + I hu, dA. external force ) (IIIB. 2) Ry oR, - J 74 54(4,5) dR internal force y - | pu. u. dR, inertial force J Ry in which F 0 on RY (IIIB.3) Substituting the incremental constitutive equation (IIIB.1) into (IIIB.2) yields: i.u uz . de. dR_ = external force [oauy, dR, + | ig di akg AUy 9 y ex (1118.4) R R - internal force y y - inertial force where Ts 6 (IIIB.5) diske ~ Sake * The°%k We observe that ds ike will possess the major symmetry if the Cia S do, which leads to a symmetric (tangent)stiffness matrix in a Galerkin/ finite element formulation ‘as will be discussed in the next section. This is due to the fact that the Truesdell rate definition of Crake in (IIIA.10) is employed. (There are other rate definitions which give a -36- symmetric tangent stiffness matrix.) Equation (IIIB.1), together with (IIIB.4) and (IIIB.5), consti- tute a very general formulation which is applicable to arbitrarily large deformations. If large rotations can be precluded, we set — * = Sijkg = CF kg 0 (I1IB.6) Consequently, ds ske = Ci ake (IIIB.7) IV. Galerkin/Finite Element Formulation and Matrix Equations IV-A. Generalized Mixed Lagrangian-Eulerian Linearized Variational Equation Combining the results that have been developed in the previous two sections, we can write the general form of the linearized mixed Lagrangian-Eulerian variational equation as follows: ; ‘ | pu, A(u, ) dR, + fotes(au,) dRy + | pu,Ac.(us) dR, R R R xX y 7 \& y y J linearized part of linearized part of convective force time derivative term + | Us 594 jke Aur ig aR, + | Us qi sky Auy 9 dR R R Nard — ) Ne ~ J linearized part of linearized part of viscosity internal force and penalty terms = | b-u. dR, + l, hu. dA, external force R oR y y ~37- - 73 5"(4,4) dR internal force y (IVA.1) - | e(u, Ju, dR, Ry inertial force in which uz. = 0 on aRy (IVA. 2) The first term corresponds to the linearized part of the time- derivative term and the second two terms correspond to the linearized part of the convective force. If we select the Lagrangian description, the latter terms are equal to zero. The third term is the tangential stiffness matrix, as described in Section III-B. The fourth term is the linearized part of the viscosity and penalty terms. We will define na d.. ijke later in this section. IV-B. Matrix Equations Background information on Galerkin/finite element formulations may be found in [18,25,26]. Review articles and current trends in finite element research can be found in [27,28,29]. The Galerkin/finite element method consists of the discretiza- tion of a region, say Rs into non-overlapping subregions called "elements." We label each element by an element number "e", e.g., the eth element is denoted by Re e=1,2,--+,Nume2. The discretization of -38- Ry also consists of Numnp "nodal points.” The position vector of the ath node, a=1,2,°°°,Numnp is denoted by Xa: The “shape function" associated with node a is denoted by No and it satisfies the relation N(x) = Sab We employ the isoparametric concept [25] that the same shape functions are used for both geometric interpolation and kinematic variables. Specifically, we assume for a typical element ("e"), with n n nodes, that: e "en e Xe = ah NaXia (IVB.1) Men e Y; = an N, Yia (IVB.2) en e Us = ay N, Ua (1VB. 3) Us = 2 N, Us, (1VB.4) us = ay No Wa (IVB.5) At "en a,e@ I 6 . . : . - ! e where Xia? Via? Usa? Usage Usa? and Usa are the a"" nodal values of t coordinates in Ry and Ry particle displacements and "mesh" displace- ments, particle velocities and "mesh" velocities, respectively. Also for programming purposes, it is convenient to introduce matrix counter- parts of the variational equations (IVA.1) and the corresponding con- stitutive equations (IIIB.1). These are summarized as follows. -39- Variational equations where n en b 2m rye cio ra dR AC y Nume2z _e! | e=1 e y T e Ba D Bh aR, Au. T Y Y “e Ba Bb Be aR, Au, T T nNibdR + | Ni haa ony BY’ + aR Tear) AR ma oy |. oN, (u ) 'y y PN, ndRy - N, 0 0 0 N 0 |= N Lo oO N, PUT 7 PHY 2 = PU] 3 PUD 7 PUD 2 = Pa 3 PUZ 7 PU3,2 43,3 (1VB.7) (1VB.8) (I1VB.9) zm ~a Hl L Na,3 a.3 -40- oOo lUlUDUlmlUlUCUCOUCUCUCDUCUO (IVB. 10) (IVB.11) (1VB.12) (IVB.13) (IVB.14) Nao2 Naat 0 Qe _ BS = 0 N33 “N32 (IVB.15) | “Na,3 0 Nal D = (6, ,]- [D; spo] (the components will be (IVB.16) ~ J defined in the next section) C147 1 Coty 2 * ©3403 N=] Cy + Coy 5 + Cally 3 (IVB.17) | 14341 * Collg.2 * 6343.3 where On denotes the mxn zero matrix, C is a 6x6 matrix of material tangent moduli whose components are given by Cry = Ci ike in which the relationships between the indices are given by: I J ; k 5 L ] 2 2 3 (IVB.18) 3 ] nN nm FP Ww DY — wOoNnNMNwWon- The ordering convention in (IVB.18) may seem strange at first, but it is quite convenient for reducing to two-dimensional and axisymmetric theories (see Part II of [16]). The initial-stress matrix, T, is defined as follows: T Te j T Te 7 r 2 6 2 6 7 9 0,8 co r!:!F 0 > j ’ Wy iy % TiN Te ts nl Z| rl a 1G rl i ; T T T | 5 it 35 T3 ! 0 Z 0 I= 5 5 0 aa. ae ee T T t T 5 6 5 5 6 ]™ 2 Bi 8 9B 2 4 a T3*Ty T2677 83 72 T= q Fi" 4 q q tT,tT, | T T T,-T att) 5 oT 1774 SYMMETRIC n : i - | m m3 Te _ 18 rn Z q T3tTy ; T%) q 4 T4tty LL 4 ~ (IVB.19) where t= T4 (IVB. 20) The D matrix is arranged to be compatible with the following ordering of strain and rotation components: 112211 22% 20°¥3322Y 23 22731 224) 222673 92434 (1VB.21) in which Vig = (Aug 5 + U5)? and (1VB.22) -43. Incremental constitutive equation The following definitions are used to derive the tangential initial stress matrix, but an improved representation can be used to integrate the constitutive equation (see Chapter 3, Sect. III). At = Ty + Sw Incremental Constitutive Equation (IVB.23) AT = {ary} (IVB. 24) € =C+c* (IVB.25) C = [C, 4 (IVB. 26) T1 4T T T 0 tf ol 2 F Truesdel I— Rate Stress - 1. 3 T2731 35 0 Matrix * = |] eee eee ee ee ee ee we eee ee eee ee ew c | (1VB.28) “Ty 0 “Ty ! Ty Te. T¢ T t TatT T 6 3.4 2 “tT eT T ' T TptT 5 2 4™°1 a ee ee F buy | Au) 9+ Aus 4 n Au en y= 2,2 = J BYA ue (1VB.29) ~ Au3 3 a=] ~™ ~ Au + Au Incremental] 2.3 32 Strain Vector Au 4 + Au, 3 ~ T) 0 “Tg = T3-T 1% ; T5 2 2 2 “To, Ts 0 S = | naeenen----2--2022- (1vB.30) 0 -T. T% 7% Tatg 72 2 2 Tp _ 72 Ty] - Ty | 2 2 Z| a ~ Auy 2 ~ Auo 4 Non - _ @ e we | Aus 3 - Au3 » = an Bo Au, (IVB.31) Au3 ~ AU) 3 | Incremental Rotation Vector ro Aus] @ e Element Nodal Displacement Aus = Aus Increment Vector (IVB.32) e L Au.3 | IV-C. Definition of di sky for Incompressible Fluid and Slightly Compressible Fluid We restrict ourselves to incompressible or slightly compressible fluids. In this way S, C*, and T can be set to zero, so that: disk = Sajke ~ Tike (1VC.1) Recall from equation (IIA.3) that .. = <-pd.. + pe . Tj PS; 2u W045) (1VC.2) For the penalty function formulation 115], we Tet -45- p= -A(u; ;) (1VC.3) where 4 > 0 is the penalty parameter. Substituting (IVC.3) into (IVC.2) yields and therefore Dry = [ds aye (Ivc.5) and f A+2u OO d r 0 07 0 u 0 0 0 0 ~ d 0 A A+2u—OO 0 0 0 0 0 u 0 L 0 0 0 0 0 Ta (Note that the § is formally the same as the matrix for classical linear isotropic elasticity.) If we consider the slightly compressible formulation, recall equation (IIA.11), t 0 + 2u WG §) (IVC.7) A time-discrete form of the pressure part of the constitutive equation is -Ap = AtpB Uk (IVC.8) from which it follows that cab) At pB + 2p Atos AtoB 0 0 0 0 0 -46- AtpB Atop 0 0 AtpB+2u Atop AtpB AtpB + 2u 0 0 ) 0 (1vC.9) 10. -47- References A. C. Eringen, Continuum Physics, Vols. I and II (Academic Press, New York and London, 1974). K. S. Pister, "Continuum mechanics--a brief review,” Lectures on Finite Element Method in Continuum Mechanics (J. T. Oden and E. R. de Arantese Olivera, eds.) (University of Alabama, Huntsville, 1973). Y. C. Fung, Foundations of Solid Mechanics (Prentice-Hall, Engle- wood Cliffs, New Jersey, 1965). J. T. Oden, Finite Element of Nonlinear Continua (McGraw-Hill, New York, 1972). J. K. Knowles, £. Sternberg, "Notes on mathematical theory of elas- ticity, Applied Mechanics (AM) 135 and 136", California Institute of Technology. L. E. Malvern, Introduction to the Mechanics of a Continuous Medium, (Prentice Hall, 1969). A. A. Amsden and C. W. Hirt,"YAQUI: An arbitrary Lagrangian-Eulerian computer program for fluid flow at all speeds," Report LA-5100, Los Alamos Scientific Laboratory, Los Alamos, New Mexico, March 1973. J. Donea, P. Fasoli-Stella, and S. Giuliani, "Lagrangian and Eulerian finite element techniques for transient fluid-structure interaction problems," Trans. Fourth Internat. Conf. on Structural Mechanics in Reactor Technology, San Francisco, August 1977. T. Belytschko and J. M. Kennedy, “Computer models for subassembly simulation," preprint. T.d.R. Hughes, W. K. Liu, and T. K. Zimmermann, "Lagrangian, Eulerian finite element formulation for incompressible viscous flows," Presented at the U.S.-Japan Conf. on Interdisciplinary Finite-Element Analysis, Cornell University, Aug. 7-11, 1978. I. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. -48- J.E.A. John and W. Haberman, Introduction to Fluid Mechanics, (Prentice-Hall, Inc., 1971). J. Chow,"Notes on Continuum Mechanics (Mate 316)," University of Iilinois at Chicago Circle. T.J.R. Hughes and J. M. Winget, "Finite rotation effects in numer- ical integration of rate constitutive equations arising in large- deformation analysis." Internat: J. Numerical Math. in Eng., in press J. O. Hallquist, "NIKE2D: An implicit, finite-deformation, finite- element code for analyzing the static and dynamic response of two- dimensional solids," Lawrence Livermore Laboratory, Report UCRL-52678, University of California, Livermore, March 3, 1979. T.J.R. Hughes, W. K. Liu, and A. Brooks, "Review of finite element analysis of incompressible viscous flows by the penalty function formulation," J. Computational Phys. 30, 1-60 (1979). T.J.R. Hughes and W. K. Liu, "Nonlinear finite element analysis of shells: Part I - Three-dimensional shells; Part II - Two-dimensional shells." Computer. Meth. in Appl. Mech. andEng., in press. T.J.R. Hughes, "On consistently derived tangent stiffness matrices, in preparation. T.J.R. Hughes, "AE/AM108: Finite element method, class notes,” California Institute of Technology. R. M. McMeeking and J. R. Rice, "Finite element formulations for problems of large elastic-plastic deformation,” Internat. J. Solids and Structures 1], 601-616 (1975). T.J.R. Hughes and K. S. Pister, "Consistent linearization in mech- anics of solids," Computers and Structures 8, 391-397 (1978). R, D. Krieg and S. W. Key, “Implementation of a time independent plasticity theory into structural computer programs," Constitutive Equations in Viscoplasticity: Computational and Engineering Aspects, AMD 20 (ASME, New York, 1976). 22. 23. 24. 25. 26. 27. 28. 29. -49- C. Truesdell and W. Noll, The Nonlinear Field Theories of Mechanics, Vol. III/3 in Encyclopedia of Physics (S. Flligge, ed.), 6pringer-Verlag, Berlin-Heidelberg-New York, 1965). M. A. Biot, Mechanics of Incremental Deformation (John Wiley and Sons, Inc., 1965). T.J.R. Hughes and J. E. Marsden, The Mathematica] Foundations of Continuum Mechanics; in preparation. 0. C. Zienkiewicz, The Finite Element Method, 3d edition (McGraw- Hill, New York, 1977). R. H. Gallagher, Finite Element Analysis Fundamentals (Prentice- Hall, Inc, 1975). G. Horrigmoe, Nonlinear Finite Element Models in Solid Mechanics (The Norwegian Institute of Technology, University of Trondheim, Norway, 1976). T.J.R. Hughes, "Current trends in finite element research," Feature article to appear in Applied Mechanics Reviews. T.J.R. Hughes, "Recent developments in computer methods for struc- tural analysis," Nucl. Eng. and Design 57, No. 2, May 1980. ~50- Chapter 3 NUMERICAL SOLUTION TECHNIQUES I. Numerical Integration of the Variational Equations The exact integration of the expressions appearing in the varia- tional equation (IVB.7 of Chapter 2) is, for all practical purposes, impossible, especially for complex, distorted elements. So instead, we employ numerical integration formulas te carry it out. For purposes of finite element analysis it is only necessary to use "sufficiently accu- rate” rules. With numerical integration we replace the integrals by a weighted sum of values of the integrand. However, the cost of numerical integration (i.e., computer time) can be quite significant. It is therefore of interest to determine the minimum integration requirement permitting "convergence", and yet to preserve the "rate of convergence" which would result if “exact integration" were used. There are many numerical integration rules which may be used. In finite element work, the Gauss numerical integration or Gauss-Legendre quadrature, as it is now frequently called [1,2,3], is well suited for Cartesian produce subdomains (e.g., quadrilaterals and bricks). I-A. Selective Integration Historically, kinematical constraints have proven difficult to deal with by the standard "displacement" finite element method. Examples of these kinematical constraints are: 1} "Nearly-incompressible" solids, in which the volumetric deformation is much much less than the deviatoric deformations. 2) “Incompressible” fluids, in which the divergence of the velocity -5]- field is zero. 3) Beams, plates, and shells based upon theories which account for transverse shear deformations in which the transverse shear defor- mations approach zero in the thin beam (plate, shell, respectively) limit. Simple, low-order elements have behaved particularly poorly. "Mixed" and "hybrid" [2,4-7,11] elements have been proposed as alternatives. The whole subject was poorly understood throughout the 1960's. In the early 1970's reduced integration procedures began to be used. However, they were also poorly understood and not generally adopted. In 1974, Fried [8], Nagtegaal et al. [9], and Malkus [10] performed important investigations focusing on the rank condition, constraint counts, and equivalence theorems, respectively. These workers created renewed in- terest in reduced/selective integration techniques and allied topics and, as a result, successful methods for kinematically constrained media could be developed within the displacement method. Within recent years, a greater understanding of this subject has been achieved. Two techniques have proved very popular. They are the reduced/selective integration procedure as advocated by Malkus and Hughes [7], and the mean-dilatation formulation proposed by Nagtegaal et al. [9]. Both techniques have proven successful, but both have their limitations. The reduced/selective integration technique is a fairly general concept and has been applied to a wide variety of finite elements. For details see [2,7,12-17]. In [12], Hughes and coworkers applied this technique to finite element analysis of plates, employing the basic -52- isoparametric elements [2]: the four-node bilinear quadrilateral and the nine-node biquadratic quadrilateral. In [13], Hughes and Cohen proposed the "heterosis" finite element for plate bending in which they used nine- node biquadratic shape functions for rotation fields and eight-node "serendipity" shape functions for displacement field, together with the reduced/selective integration technique. In [14] Hughes and Liu suc- cessfully extended the technique to nonlinear finite element analysis of shells as will be discussed in more detail in the next chapter. In [15,16] Hughes and coworkers applied the technique, together with the penalty function formulation, to finite element analysis of incompressible viscous flows and Lagrangian-Eulerian descriptions for incompressible fluids. They will be discussed in detail in Chapter 5. In [17] Hughes and Prevost employed the same ideas in nonlinear quasi- static finite element analysis of soil via the computer program. DIRT II. The mean-dilatation formulation is a more special technique which has been used mostly in the context of basic isoparametric ele- ments. The two approaches are identical for the four-node quadrilateral in two-dimensional plane strain and the eight-node brick in three-dimen- sional analysis. However, the mean-dilatation formulation appears to behave somewhat better in axisymmetric analysis than does the correspond- ing reduced/selective integration procedure [9]. The rationale behind these two techniques has been discussed at length in [2,7-10,12-17], so we only provide a brief summary here. The reduced/selective integration procedure for the continuum elements case goes as follows: The "stiff" term (dilatation term) may be segregated from the remainder (deviatoric term). A reduced quadrature rule is then used on the stiff term to lessen the constraint -53- ("locking") and a normal, quadrature is used on the remainder such that correct rank is maintained. If the segregated terms are tensor invari- ants, that is, the resulting matrices transform properly with respect to rotation of the global reference frame [12], then the resulting arrays are similarly invariant. An example of selective Gauss-Legendre inte- gration rules for two-dimensional isoparametric Lagrange elements is presented in Figure IA-1. Extension to three-dimensional isoparametric Lagrange elements is straightforward. One drawback of the reduced/ selective integration procedure is that triangular elements and "seren- dipity" quadrilaterals generally exhibit inferior behavior [18]. In the context of the general material laws, such as equation (IVB.23) of Chapter 2, explicit segregation of the contributions to the matrix equations into dilatation, deviatoric, and coupling terms proves very inconvenient. In principle, a decomposition of this kind is pos- sible, however, it is not at all clear how to treat the coupling terms. Also the element implementation is awkward. Clearly an alternative general scheme is called for. Hughes proposed a generalization of selective integration pro- cedures to anisotropic and nonlinear media [19]. Instead of separating the d. (in eq. (IIIB.5) of Chapter 2) into dilatation, deviatoric, TjkR and coupling terms, the matrix BY appearing in equation (IVB.7) of Chapter 2 is modified. A summary of the procedure is as follows; for details the reader is urged to consult [19]. The original BY is dil a divided into B and poev which are defined as: dil d BY =B + Bory (1A.1) ~ ~ -54- ‘sjusUlajTa eBuersey] s11y9uUTeredost [euotsuaulp “2 IO} sent uoTeIZojuy arpuasey-ssneH aayoejag *|-y{ bly vx exe Zxz Wwia}-77 GX cxzZ jyuiod | W419}-Y IIGNIIG d1,D4ponbiq IDAUIIG SUO!JOUNS adpys -55- where - + B, = Bn Bg 0 0 0 B B B . 7} 2 3 dil _ 7 |] 0. pf a1 (1A.2) B, Bo By 0 0 0 L 0 0 0 J B, 0 0 B, Bo B3 Bo By 0 0 0 0 dey 0 Bo 0 By By B. Boo | woccttecceercccee Tgp ocrccrccceeccrecsee (IA.3) 0 0 B. B. Bo B, 0 B. Bo 0 0 0 LB 0 8 | L 0 0 0. We construct a new BY called Br (which replaces BY) by: BY = pdev 4 pall : (1.4) ~a 4a ~a in which Bo"! is the "improved" dilatational contribution and is defined Py: f BBR BST) 0 0 0 ee ac (1A.5) By By By 0 0 0 - 0 0 0 + The resulting BY is - 4 B, Bg SB B, By 0 : B, By By BY Sf] ee ee eee eee (IA.6) B, Be Bg 0 By By L By 0 By where _ ~ By = (B, - By) / 3 B, = By + By B. = (B, - B,)/3 6 2 2 (1A.7) By = By + Be Be = (B, - B,)/ 3 By = By + By J With By defined, the stiffness matrix and internal force vector are integrated within one “poLoop’ [20], using a single integration rule. pai! tegration of all orders; (2) mean-dilatation and higher-order general- can be constructed to be equivalent to: (1) selective in- izations. Similar procedures have been employed by Hughes and Liu [14] in nonlinear finite element analysis of plates and shells (see also Chapter 4). As an example, let us pick Be such that it is equivalent to selective integration. A quadrature rule is specified to integrate the element stiffness matrix and internal force vector. This rule is con- sidered the "normal" one. Introduce a “reduced” rule; let Neng and Ey be the numbers of integration points and the corresponding locations -57- for the "reduced" rule. A specified set of shape functions, N,'s, is de- fined with the nodal points equal to E, s(i-e., NCE) = 6.5). The general form of B is int B,(é) = } wi ky a N,(E) Bap (1A.8) The equivalent of selective integration is achieved if Bay = B.(E,). The specialization of equation (IA.8) to mean-dilatation formula- tion and other generalizations has been discussed by Hughes [19], and will not be discussed here. A flow chart of the computer programming aspects is presented in Table IA-1. Table IA-1. Flow Chart of Generalizations of Selective Integration Procedures 1. Loop over each element. 2. Pick up position vector arrays. 3. Define the shape functions Ns (E) and derivatives B. o(E) for l<k < Ne igs and store. 4. Loop on the “normal” integration points (E,)- 5. Compute the "normal" shape functions, N.(Eo)s and corresponding derivatives B.(E,). 6. Define B.(E,) using equation (IA.8). 7. Define BYE.) using equations (IA.6) and (IA.7). 8. Form stiffness matrix and internal forces. 9. If &< Neo GO TO 4 t If £>N, 60 T01 t -58- I-B. "Upwind" Finite Elements by Modifying the Standard Gauss-Legendre Rules As can be seen from equation (IVB.7) of Chapter 2, the second and third terms constitute a nonsymmetric operator. The Galerkin finite element method has been successful in structural and solid mechanics in which the operator is a symmetric one. Until recently, nonself- adjoint cases (nonsymmetric operators) have not been well understood. If we use Gauss quadrature on the convective terms appearing in the above-mentioned equation, spurious spatial oscillations are exhibited for some high Reynolds number flows. These noisy results can only be corrected by mesh refinements. However, in most problems of engineer- ing interest, convection is dominant and the geometry is so complex that severe mesh refinement is not economically feasible. To circum- vent these problems, "upwind" finite elements [15,16,25-29] or Petroy- Galerkin methods [22-25] are used as alternatives to generate stable finite element approximations to highly convective flows. As observed in the finite difference literature, "upwind" dif- ference schemes [21] may be employed to preclude such oscillations, but due to excessive numerical diffusion, the accuracy of the solution may degrade considerably. Among the proposed upwind-type finite element procedures , some schemes are subject to the same criticisms. However, some have been shown to be superior (see e.g., [28,39,40]). In fact, for some convection-dominated flows, the basic Galerkin/ finite element method is adequate, and a great deal of success and progress has been reported in the finite element literature (see e.g., [30-36]). Nevertheless, based on the "Comments" and "Reply to the Comments" made by Gartling [37] and Zienkiewicz [38], respectively, and -59- also from the numerical example shown in [15], basic Galerkin/finite ele- ment methods exhibit spurious behavior wherever there is a downstream essential, or "hard" boundary condition, in convection dominated situa- tions. The first method used is based upon numerical integration rules which are simple, efficient, very easy to implement, and yet effective. Instead of using the standard Gauss-Legendre rules, a modified "upwind" integration rule was employed. The positions of the integration points are determined from the “element Reynolds number." Since it has been adequately documented elsewhere [see 16, and Sec. 4.3 of 15], this modi- fied "upwind" integration rule will not be discussed here. However, we will show some numerical examples to demonstrate the effectiveness of the proposed theories in Chapter 5. The performance of the "upwind" scheme used was found to be sat- isfactory for a wide class of slow problems, as will be discussed later. However, for some time-dependent convection-dominated flow problems it was found that the "spurious crosswind diffusion" produced by the upwind scheme degrades the accuracy of the solution [see 28,39]. A multidimen- sional upwind scheme with no "crosswind diffusion" was designed and shown to be superior to the quadrature scheme [28,39]. Computational fluid dynamics is still a relatively young science so it is still too early to evaluate the relative merits of the various proposed schemes. Although these schemes differ considerably in concept and implementation, they are "designed" to achieve similar ends. There- fore, it is beyond the scope of this thesis to design the "ultimate" scheme for computational fluid dynamics. The interested reader may profitably consult [25-39] and references cited therein for further in- formation on this topic. -60- While the objective of this thesis is development of finite element procedures for fluid-structure interaction, it is worthwhile to implement new ideas and examine relative merits. Therefore, the upwind scheme with no "crosswind diffusion" proposed in [39] will be used for the examples shown in Chapter 6. However, for fluid-structure interac- tion problems in which the mixed Lagrangian-Eulerian formulation is employed, the scheme in [39], which is restricted to the Eulerian case, has to be modified. Following [39], we modify the equation of motion in the mixed Lagrangian-Eulerian description (equation IB.10) of Chapter 2) to: o(uy) + pes(uy) 5 = (ty5 * 45) 5 (1B.1) where Dj is an artificial viscous force and is determined by: Ds; = Usk Way (1B.2) in which G, = ¢&/ llcll (1B.4) and u is defined as follows: Ces ) ( ) Ll =. E. cp: h,.)/N IB.5 gf) 71 Ei 1 SD where E. ; coth Oe go Vor; or approximately a a, 4/3 “35%,5 53 5 (1B.6) SQN Op, [oe | > 3 in which Oey = PCr; he, /2n (1B.7) and _ Coy = Oey * c (1B.8) The meaning of ne and er; is illustrated in Fig. IB-1 for the four-node bilinear quadrilateral. The corresponding variational equation for (IB.1) can be written as follows: fo Cou, (ay) + ays 5] aR, R y + | p(u, +4; ec, Us dR, y = | b.u. dR + I, hau. dA (IB.9) R y aRy in which — . g u, = 0 on aR (1B.10) and u, = — (1B.11) i plici] “isk Sk In [28], it is suggested to use "consistent weighting." For the intended applications it was thought to be sufficient at this stage of development to neglect consistent weighting. Nevertheless, it is anti- cipated that a consistent weighted formulation (Petrov-Galerkin scheme) would lead to improved results in transient analysis. (This is being -62- —e Xx) Fig. IB-1 -6§3- undertaken in the thesis of Brooks [40].) I-C. Mass Matrix If we use the standard Gauss quadrature on the first term of equation (IVb.7) of Chapter 2, the resulting matrix is called the con- sistent mass matrix. It is symmetric and positive-definite. However, for several computational reasons, "Jumped" mass matrices are desired. For two- and three-dimensional, rectilinear cases we employ product Lobatto integration rules [2] on Lagrange elements which gives rise to a diagonalized mass matrix. The first two Lobatto rules are the trapezoidal and Simpson's rules, whose products are appropriate for the 4- and 9-node elements, respectively, in the 2-dimensional case, and whose triple products are appropriate for the 8- and 27-node elements, respectively, for the 3-dimensional case. In the axisymmetric case, zero masses along. the z-axis, due to the factor r (radius) appearing in the volume element, would result if Lobatto integration rules were employed. This can cause difficulties in transient analysis. Thus a "row sum" technique is employed to cir- cumvent the problem. The consistent mass matrix term is e _ T Mob | pn. Ne dR, (1C.1) Ry The row sum matrix is defined to be N e _ ren e - ma Mb no sumona, a= 1, -++sNQ, (IC.2) and thus we arrive at the row sum mass definition: -64- e Mb = Sab [. en, dR no sum on a, a Tyree sNa, (IC.3) R YN using the fact that J" N, = 1. We employ Gauss quadrature to evalu- a=] ate the integral in (IC.3). For direct integration procedures, the type of mass matrix that should be used depends on the method of temporal integration as will be discussed later in this chapter. It is interesting to point out that consistent mass matrices tend to overestimate the frequencies, whereas lumped masses underestimate the frequencies. However, mass lumping is still a controversial issue in fluid mechanics, due primarily to the results of Gresho et al. [41]. Il. Implementation of Lagrangian-Eulerian Finite Elements Since all the volume and surface integrals in the variational equa- tions and thus the matrix equations are in terms of Ry and OR» respectively, it is important to specify what Ry is in our finite ele- ment discretization. We wish to think of the region Ry as our finite element mesh in the referential domain. Ry is the image of RY at time t é [t,t]. Hence, $ is the "motion" of the mesh. In general d is arbitrary. The selection of $ is generally made depending on the use- fulness for the intended application. For our purpose, fluid-structure interaction, it is useful to specify u's or equivalently to specify u and thus 6. So u is the mesh displacement, u' is the mesh velocity, and u" is the mesh acceleration. If we select u’ = 0, we have the Eulerian description; If we select @' = u, we have the Lagrangian description. If we select a -65- between 0 and u, we have a mixed description which is employed in Chapter 6 for fluid-structure interaction problems. We will employ the idea of split descriptions and relate u' to u 1 > @o> @a> and Uy following classical ideas [42]. In this way the unknowns u, and u" can be solved by temporal numerical integration [2,3,46], and $ is then uniquely specified. For further generalizations see [43-45]. ~ II-A. Generalized Description Using Isoparametric Finite Elements In this section it is assumed that unless otherwise indicated, there is no sum on repeated i indices. Let us introduce a "Lagrange- Euler parameter array” defined by: a= [a,j] = [oi 555 5] (IIA.1) where a, € [0,1], for i=1,--+.Ncp. We relate u' to a and u by: x, d ‘sd ; a(u; +x.) j=l k=1 k . #3 E or yn fixed where Cy. = Ups Uy = (5; - 0,5 14, (1IA.3) We assume that the a5. 'S are constants in the present developments. As an example, consider a two-dimensional mesh, and let 11 = 0, Gog = 1, {1IA.2) becomes uy = 0 and . . a(ti5 + Xp) (IIA.4) Us = Uy - U ef 1 8% be on n Fixed -66- If a.. = O we have u' = Q, the Eulerian description. Ifa. = 6.. Ww ij “ ~ ij ij e have u' = Us the Lagrangian description. We may have split or mixed Lagrangian~Eulerian descriptions as in IIA.4. If Os takes on values strictly between 0 and 1, we have a very general description which is suitable for bubble dynamics and other related applications. See [43-45,47,48] for further details. However, for our applications we restrict ourselves to the first three cases. For finite element implementation, we assign a value of @ to each nodal point. Further, we employ isoparametric interpolation and let Nen — e 7 = oe Oe = , Nia i=l, Neg (IIA.5) a=] However, in attempting to evaluate the derivatives in (1IA.2), by virtue of the fact that we are using c°-continuity shape functions, dis- continuities will result across element interfaces. Consequently, a “weakened” interpretation of (IIA.2) is necessitated. To this end we may introduce Galerkin approximations of (IIA.2). Let v= {va qh (1IA.6) v= {vig} (IIA.7) A= [a 9] (11A.8) where Vs v, and A are the generalized global vectors of nodal mesh velocities (u'), material velocities (u), and Lagrange-Euler parameters -67- (a diagonal matrix), respectively. The subscript A takes the values 1,2,-:+,Numnp. The subscript i takes the values Taree Nae With these definitions, the Galerkin equations may be written as: f= aly LS) (118.9) where Nume& L= J) f nTy dR (IIA.10) e=] Re X Nume2 Non n a) [x At, Nae, (aC, + 8°,) dR, (11A.11) ~ e=] Re a=1 3 x Here Re is the referential region of the eth element. We employ Lobatto rules to evaluate the integrals in (IIA.10) and (1IA.11). The resulting matrix L is symmetric, positive definite, and diagonal. One drawback of using Lobatto rules is that only Lagrange elements can be used, since "Zeroes" or "negative" weighting will result if other element types are employed. In the axisymmetric case, we employ area average instead of volume average to avoid zeroes along the z-axis. If Gauss quadra- ture rules are employed to evaluate (IIA.10) and (1IA.11), L is non- diagonal though banded; also, a"diffusior' of the description (Lagrangian- Eulerian) may take place [16]. III. Numerical Integration of Rate Constitutive Equations The incremental constitutive equations in Chapter 2 (IVB.23) is a linearized expression. In a numerical formulation, the stresses wil] be calculated at all Gauss quadrature points which correspond to a set of material particles. As we have noticed before, the first term on -68- the right-hand side of the mentioned equation represents the material response due to deformations, whereas the second term accounts for rota- tional effects, which are uniquely specified by "objectivity." However, standard time-discretization procedures [49,51,57,58,60] when applied to this incremental constitutive equation typically only achieve objectiv- ity in the limit of vanishingly small time steps, which is not economical in practical computing. If one employs the polar decomposition to separate deformation and rotation effects [59], much larger time steps can be accommodated. However, it is too time consuming to extract the eigenvalues and eigen- vectors [52] at each Gauss point. Hatiquist [50] has shown, in the context of two-dimensional analy- sis, how to maintain objectivity for larger time steps. However, there are some awkward aspects of the implementation, and the generalization to three dimensions is not apparent. Hughes [54] proposed an improved algo- rithm for integrating rate constitutive equations in large-deformation analysis. The algorithm is fairly simple and is shown to be objective with respect to large rotation increments. Therefore we employed the mentioned technique in [54] to account for the rotational term (i.e., Sw) and the various standard numerical integration algorithms [53,55,56] to account for the material response part (i.e., Cy). III-A. Numerical Algorithm for Rotational Term We represent the position vector of a particle at time t,= [t,t] by y", ~ Let At be the time increment. The position vector at -69- ] t,t At is denoted by y"* which may be written as a function of y" and the time step; hence, + y" 7] n+l ( = yl (y" at) (IIIA.1) The particle displacement increment over the step is gs = uetl yn (IIIA. 2) ~ Here yl and u" are the approximated values of u at time t, + At and time tye respectively. Consider the following one-parameter family of configurations: n y to (1 -a)y" $ ay" (IIIA. 3) in which a-€ [0,1]. Let G denote the gradient of 6 with respect to yr, In component form _ nto, G.. = 365 /3Y5 ij (IIIA.4) The incremental strain vector .y, and the incremental rotation vector, w, may be defined in terms of G as follows: y= (646!) /2 (IIIA.5) fl * 363 (6-6) /2 (IIIA.6) The numerical algorithm for integrating the constitutive equation is as follows: vl 2 OTs at (IIIA.7) mT sgt! (IIIA.8) -70- Q-1+ -aw) yo (IIIA.9) At = Cy (IIIA.10) The expression for Q is obtained by applying the generalized midpoint rule (see e.g., [61]) to the generating equation dQ/dt =wQ If a = 1/2 it is a second order accurate algorithm. For further details of the properties of the algorithm, consult [54]. III-B. Numerical Algorithms for Integration of Eq. (IIIA.10) As already emphasized, the incremental equations derived in Chapter 2 admit a wide variety of material constitutive models. However, the application of these equations to physically nonlinear problems requires detailed knowledge of the material characterization, that is, C and C, The constitutive algorithms one should use then depend solely on the definition of C. Therefore we will consider numerical algorithms for linear elasticity, the Saint-Venant-Kirchhoff nonlinear elastic model, the Krieg-Key combined isatropic-kinematic hardening plasticity model, incompressible viscous fluids via the penalty function formulation, and the slightly compressible fluid. These are sufficiently general for the developments herein. ]. Linear Elasticity The C matrix in this case is constant and is given by: -71- Tet 0 A A OCT 0 yo 60 0 0 c=] % 8 An A OO (111B.1) ~ h O r AMY O 0 0 o 0 0 uw 0 Lo 0 0 0 0 4. If the Truesdell rate is neglected, At can be computed easily by taking advantage of sparseness. Otherwise, a general matrix multiplication “routine” is employed. 2. Saint-Venant-Kirchhoff Nonlinear Elastic Model The C matrix is defined by the following algorithm; a) Compute the deformation gradient matrix F (Fi = a(u, +x) / ax5) b) Compute determinant of F; J = det(F). c) Calculate the Finger tensor B = FF! d) Define an array B(1-6) with: B(1) = B(1,1), B(2) = B(2,2), B(3) =B(3,3), B(4) = B(1,2), B(5) = B(2,3), and B(6) = B(3,1). e) Define an array A(1-21) using the following algorithms: K=1 | el. LOOP ON I= 1,6 BB = B(I) e2. LOOP ON J = I,6 A(K) = BB*B(J) K = K+] IF(J.LT.6} GO TO e2. IF(I.LT.6) GO TO el. f} Define constants C1 to C5 Cl = A/J C2 = u/d C3 = C2+C2 C4 = C1+C2 C5 = C1+C3 ~72- g) With the shorthand notation AI=A(I), the C matrix is then defined as follows: C5 «AI C5«A4 C1*A2 C1*A3 C1*A5 C5*A6 +C3*A16 | +C3*A21 | +C3*A18 C4*A16 C5*A9 C1xA13 C4xAl7 C4xA18 +C2xA2 +C3%A20 | +C2*ATT | +C2*A5 C5%A7 C1*A8 C5*A10 C1*A11 +C3*A19 +C3%A17 C = ~ C5*A12 C5«Al4 C5*A15 C4xA19 C4*A20 SYMMETRIC +C2*A8 +C2*A13 C4*A21 +C2*A3 L J (IIIB. 2) At can then be computed using a general matrix multiplication routine. 3. Krieg-Key Combined-Isotropic-Kinematic Hardening Plasticity Model The For this model the radial return method [60] is employed. numerical algorithm is as follows: a) Compute the elastic "trial" stress increment as described in the linear elasticity section. Update the elastic trial stress and call it t®, using (IIIA.7). -73- b) Obtain the deviatoric part of (1° - al): zt -a rs rs rs c) Compute the second invariant, Jos 2 2 2 = deviatoric part of (xt n) (IIIB.3) Jo= & (CET ~ Eby)” + (Gy = EGg)” + (E53 ~ Ey) ]# ppt E43 Eg (1118.4) d) If (3J, - (k")2) is less than or equal to zero, vl = t's go to h. e) Compute plastic strain increment Ad, AdP = (Z- k")/(3u +B) (IIIB. where Z = v¥3d5 f) Compute the following quantities: AS = 3y Ad?/Z (IIIB. Aa = (1-8)B Ad?/zZ (IIIB. +] n+] ; position of yield centers, a 3 ~ g) Radial return for stresses, a" and yield stress Kel ptt. tL ag ct (IIIB atl = aM + naet (I1IB KetT - gD a gpa dP (IIIB. h) If the tangent stiffness matrix is not required, exit. i) Compute the C matrix for tangent stiffness calculation as follows: 5 “ eo + ct | ou (IIIB. Ne" (3usB) 6) 7) .8) 9) 10) 11) «means "is replaced by." -74- r At2u + _ 2 “15 X-&63 | A &Sa]) bye] ~ F186 7 2 ~ &&3 | ~ &9bq - 6965 | ~ 98 At2u 2 | *~&3&q | ~ &3& | ~ &3& C = ~ A+2u 2 ~ & abs | ~ £455 SYMMETRIC 2 2 uo &6 ” a (IIIB. 12) In (IIIB.12) the shorthand notation er = big is employed (as described by (IVB18 ) in Chapter 2. 4. Incompressible Viscous Fluid via Penalty Function Formulation By virtue of the fact that the incompressible viscous fluid used is isotropic, equations (IIIA.7) through (IIIA.10) are replaced by: tte Ma at (111A.13) ~ (IIIA.14) AT ~ Cy where C is defined by (IIIB.1) and -y is defined by (IIIA.1) through (IIIA.5) in which (IIIA.2) is replaced by: -75- p= umtl_ gn The Lamé constant A appearing in (IIIB.1) is interpreted as the “penalty parameter” » which has been described in detail in [15]. Of course, due to the sparseness of C in (IIIB.1), At can easily be computed with- out employing a matrix multiplication routine. 5. Slightly Compressible Viscous Fluid We assume the fluid is slightly compressible. Equations (IIIA.13) and (IIIA.14) are employed. However, we divide the At computation into two parts, the viscous part and the pressure part. The viscous part can be computed using the procedures described above for the incom- pressible viscous fluid with A set equal to zero. The pressure part is computed as follows. We assume the pressure at step n, p', is known and stored at each Gauss point. The divergence of the material velocity field, 8,is computed at the nta geometry, that is: a = au) /a(yF"%) sum on i= 1,+-+ Nop (IITA.15) in which uM = (T= a)u + wl!» ete. (IITA. 16) +] Then the new pressure po" is: pl = 5" - at ppe (IIIA. 17) Since we approximate u(t.) by u", the divergence part of the linearized expression A div u (appearing in the 5th term of eq. (IVB.7) in Chapter 2) is approximated by A div uF = @ div Au which is consis- tent with the pressure computation in which we evaluate 6 at the nta -76- geometry. For implementational purposes, we include o in the tangent matrix (eq. (IVC.9) in Chapter 2) and _ — — ~ At2u 0 A A 0 0 0 u 0 0 0 0 . XO wu F 0 0 i _ _ (IIIA.18) ~ r 0 r At2p 0 0 0 0 0 0 uu 0 | 0 0 0 0 0 wo in which i = AtpBa. IV. Transient Analysis Techniques for the Variational Equations There are two strategies used in dealing with the transient prob- lem: (1) modal superposition, and (2) direct integration [2,3]. In modal superposition methods, the variational equations or matrix equa- tions are diagonalized by finding the eigenvalues and eigenvectors. These eigenvalues and eigenvectors correspond to the natural frequencies and modes. In direct integration the matrix equations are integrated using a numerical step-by-step procedure without transformation of the equations into a different form like modal superposition. In linear analysis, the choice of method depends on the frequency content of the load and on which portion of the frequency response is of interest. For some particular cases, a combination of both methods is very effec- tive [62]. In modal superposition procedures, the most critical and time- consuming aspect of the computation is the determination of the eigen- values. There are a number of eigenvalue methods used in practice [1,3,46]; each one of them has its own merits. Modal superposition is -77- used very little in nonlinear analysis. It appears that it is not suit- able for path-dependent material and finite deformation problems. For most nonlinear problems, direct integration is the only ef- fective method in transient analysis. Since we are dealing mostly with nonlinear problems, only direct integration techniques will be discussed. Basically, there are two general classes of direct integration methods: (1) implicit integration, and (2) explicit integration [2,3, 46]. Both methods are developed from difference formulas that relate the accelerations, velocities, and displacements. Implicit algorithms, some of which are unconditionally stable, permit large time steps, but the cost per step is high and storage requirements tend to increase dramatically with the size of the mesh. On the other hand, explicit algorithms, which are conditionally stable, tend to be inexpensive per step and require less storage than implicit algorithms; but numerical stability requires the size of the time step to be inversely propor- tional to the highest frequency of the mesh. Belytschko and Mullen [63,64], Hughes and Liu [65,66], and Park et al. [67-69] proposed implicit-explicit nodal partition, implicit- explicit mesh partition, and staggered solution procedures, respec- tively, to circumvent these difficulties. Other procedures in transi- ent analysis algorithms have also been proven to be successful [70-83]. Nevertheless, a generalization of al] these ideas in one package, though desirable, requires further research in theoretical and imple- mentation aspects. ~78- IV-A. Temporal Integration of the Variational Equations In this section we are going to employ the new family of implicit- explicit algorithms developed in [65,66] for the temporal discre- tization of the variational equations. Although there are many dif- ference formulas that can be used in time discretization [2,3], based on previous studies, the Newmark family of methods [84] appears most suitable for fluid-structure interaction problems. For discussion purposes, let ds V,3 and a be the nodal vectors of mesh displacements, mesh velocities, and mesh accelerations {i.e., Un» us u") respectively; likewise, das v,> and a, are the material dis- placements, material velocities and material referential accelerations (i.e., Uns uns uy), respectively. The solution at time step n+] is determined by the following equations (in order): i=0 (i is the iteration counter) (IVA.1) (i) 2 4 . dit] ~ qn] (IVA. 2) (i) _ Yatt ~ “ntl predictor phase for (IVA. 3) (i) material displacements, ant] = 0 \ velocities, and ac- (IVA. 4) ~ ~ celerations ~ _ 1. 2 anal = gd, * At n*& B)at on (IVA.5) Vat] 7 Yq t (ly) Atay J (IVA.6) a(i) _ % ‘ diel ~ dat] (IVA.7) a(i) _ & predictor phase for Vv = V mesh displacements, . Ata = Via h displ (1VA.8) ~ ~ and velocities x 4 A ] On diay = a, * At Vit (g- Bat ay (IVA.9) Vat] = Yn t (1- y)At ai, (IVA.10) = + lt AF* D-N ~C vas | ext Pat Pe = -79- AF* momentum equations m+ yat on (G61) vl) ali) + at Ky(ght)) + eat?K (vl) aC? gG)) ext 5 ead (l)) + Fed) ~ 1g) al) - ww Olah) pint) (4) 34) (id) ~ Fae net Gnd Sn Nume2 T ) | oN'N dR mass matrix e=1 Re ~™ Y y Nume& ) | (8')" pelce!) + £4(B2)] aR e=1 pe ~ * ~ * y ¥ X J tangential convective matrix Nume2 TA ) | BY' D BY dR pressure and viscous e=] Re ~~ y tangent matrix y Nume 2 T = } | B'DB dR tangent stiffness matrix e=1 pe ~~~ y discrete nodal applied forces Nume & T )} | Nb aR body forces el te~~ (IVA.11) (IVA.12) (IVA.13) (IVA.14) (IVA.15) (IVA. 16) (IVA.17) (IVA. 18) (IVA.19) -80- Nume& T Fe = | Nvh dA surface forces (IVA.20) ~ e=] gre ~~ oY y Nume 2 T (4) I = Jj pN ane] dRy inertia forces (IVA.21) e=1 ne y Nume& T No = ; | pN'n aR, convective forces (IVA. 22) e= e RY : Nume& ‘ int _ y! (i) ; Fat] L [ BY ta] dR, internal forces (IVA. 23) RY a(i#1) = ali) + na ) (1VA.24) ~ ~ ~ corrector phase (i#1) _ ~ (i+1) for material Vv = Vv + yAt a displacements, (IVA. 25) “nt “ntl ~nl velocities, and | | accelerations it]) _ x 2_(itl dnt] ~ qt + BAT an 44 y, (IVA. 26) a(it]) _ (i+1) -1 mesh velocity Vat Ana] -t S) representation (IVA.27) A(it]) _ ~(it1) * att > Viney “Vneq) / At corrector phase (IVA. 28) ( ’ for mesh a{itl) _ 4 2 alitl accelerations and ae diay * BOt™ any displacements (IVA. 29) The notations used in equations (IVA.74) through (IVA.23) have been defined in Chapter 2 (equation IVB.7). The numerical integration for if) that appears in equation (IVA.23) has been discussed in the pre- vious section. The mesh velocity update formula, equation (IVA. 27), has been discussed in Section II of this chapter. -8]- Let I be the total number of iterations. If I=0, the algorithm is a predictor-corrector algorithm. If more iterations are required (i.e., I>0), it is a predictor-multicorrector algorithm. The generalizations of (IVA.1) through (IVA.29) to implicit, explicit, implicit-explicit mesh par- titions, implicit-explicit operator splitting, dynamic relaxation, and quasi-Newton updates for symmetric and nonsymmetric equation systems have been discussed in detail in [66,72,79]. The stability, accuracy, and convergence have been discussed in [65,72,73]. For detailed analyses of the above algorithms, the reader should consult [75,66,72,73,79] and references cited therein. IV-B. Computer Implementation of Equations (IVA.1) through (IVA. 29) We assume the initial data doe Va? 452 da. Vo? ag are all given. If the mixed description is used, the global array L is computed accord- ing to equation (1IA.10) and stored. For computer implementation pur- poses, let "I" and "E" denote “implicit" and "explicit" element groups. If the entire mesh is treated explicitly, E M* = M which is a diagonal matrix (IVB.1) ~ apt = (EYE + Be + eg = (BME (18.2) The terms onthe right-hand side of equation (IVB.2) are defined by equation (IVA.18) to (IVA.23) with i set to zero. If the entire mesh is treated implicitly, M™ is defined as in equation (IVA.12) and equa- tions (IVA.13) through (IVA.17). M* can be reformed and factorized whenever necessary. However, (aF*)! is recomputed for every iteration. If a linear implicit/nonlinear explicit splitting is desired, M* is -82- defined as MY = M+ yat KI + pat’K’ (1vB.3) wk. : * where M is constant. and 1S computed once and stored. The corre- sponding AF™ is recomputed according to equation (IVA.13) for every I E (pexty! + (pextye - pext iteration. We notice that M + M” = M, and With these definitions in mind, the solution procedures are as follows: a) Loop on number of time steps, N. b) Define predictor 'material" values according to equations (IVA.1) through (IVA.6). c) Define predictor "mesh" values according to equations (IVA.7) through (IVA.10),if required. jE d) Form pext and (AF*)~ as described previously and store. ~nt] e) Loop on number of iterations, i=1,---,I. f) Reform and factorize M* if required, where M* = cm)? + (M*) . g) Form (aF*) | and AF™ = (aF*)? + (ar*)®. h} Solve for Aa according to equation (IVA.11). i) Update material displacements, velocities, and accelerations ac- cording to equations (IVA.24) through (IVA.27). 3) Form S$ according to equation (IIA.11) if required. k) Update mesh velocities according to equation (IVA.27), if required. 1) Update mesh accelerations and mesh displacements according to equa- tions (IVA.28) and (IVA.29), if required. m) Check convergence,if required. If convergence check is required then for ||Aa|| <e go to o. n) If convergence check is required and i=I stop. If i <I, go toe. . (i+1) 0) t, 1s replaced by Thal p) q) ~83- qgfitt) yt) | and alt 1), respectively. d,s int 2 a, are replaced by doi. “s Viyy a i+1) 9(itt) 4 a( dvs Vas a, are replaced by d ntl? “ntl? and If N < maximum number of time steps, go to a. a 1) antl > respectively. 10. 1. 12. -84- References C-E. Froéberg, Introduction to Numerical Analysis, 2nd edition. (Addison-Wesley Publishing Co., 1969). 0. C. Zienkiewicz, The Finite Element Method, 3rd edition. (McGraw- Hill, 1977). K-J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis. (Prentice-Hall, Inc., 1976). G. Horrigmoe, “Nonlinear Finite Element Models in Solid Mechanics" Report No. 76-2, The University of Trondheim, Norway, 1976. G. Horrigmoe, “Finite Element Instability Analysis of Free-Form Shells," Report No. 77-2, The University of Trondheim, Norway, 1977. T.H.H. Pian, "Derivation of element stiffness matrices by assumed stress distributions." AIAA Jour. 2, 1333-1336 (1964). D. S. Malkus and T.J.R. Hughes, "Mixed-finite element methods-- reduced and selective integration techniques: A unification of con- cepts," Computer Meth. in Appl. Mech. and Eng. 15, 63-81 (1978). I. Fried, "Finite element analysis of incompressible material by residual energy balancing," Internat. J. Solids and Structures 10, 993-1002 (1974). J. C. Nagtegaal, D. M. Parks, and J. R. Rice, "On numerically ac- curate finite element solutions in the fully plastic range," Computer Meth. in Appl. Mech. and Eng. 4, 153-178 (1974). D. S. Malkus, "Finite element analysis of incompressible solids," Ph.D. thesis, Boston University, Boston, 1975. T.H.H. Pian, "Hybrid models," in Numerical and Computer Methods for Solid Continua, $. J. Fenves et al. (eds.) (Academic Press, New York, 1973). T.J.R. Hughes, M. Cohen, and M. Haroun, "Reduced and selective integration techniques in the finite element analysis of plates," Nucl. Eng. and Design 46, 203-222 (1978). 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. -85- T.J.R. Hughes and M. Cohen, "The 'heterosis' finite element for plate bending," Computers and Structures 9, 445-450 (1978). T.J.R. Hughes and W. K. Liu, “Nonlinear finite element analysis of shells: Part I - Three-dimensional shells; Part II - Two-dimensional shells", Computer Meth. in Appl. Mech. and Eng., in press. T.J.R. Hughes, W. K. Liu, and A. Brooks, "Review of finite element analysis of incompressible viscous flows by the penalty function formulation," J. Computational Phys. 30, 1-60 (1979). T.J.R. Hughes, W. K. Liu, and T. Zimmermann, “Lagrangian-Eulerian finite element formulation for incompressible viscous flows," U.S.-Japan Seminar on Interdisciplinary Finite Element Analysis, Cornell University, Ithaca, New York, August 7-11, 1979. T.J.R. Hughes and J. H. Prevost, DIRT II - A nonlinear quasistatic finite element analysis program,"California Institute of Technology, 1979. D. S. Malkus, "A finite element displacement model valid for any value of the compressibility," Internat. J. Solids and Structures 12, 731-738 (1976). T.J.R. Hughes, “Generalization of selective integration procedures to anisotropic and nonlinear media," Internat. J. Numerical Meth. in Eng., to appear. P. Cress, P. Dirksen, and J. Graham, Fortran IV with WATFOR and WATFIV, (Prentice-Hal], Inc., 1970). P. J. Roache, Computational Fluid Dynamics (Hermosa Publishers, Albuquerque, New Mexico, 1976). I. Christie, D. F. Griffiths, A. R. Mitchell, and 0. C. Zienkiewicz, "Finite element methods for second order differential equations with significant first derivatives," Internat. J. Numerical Meth. in Eng. 10, 1389-1396 (1976). D. F. Griffiths, "On the approximation of convection problems “n fluid dynamics," Research paper No. 317, University of Calgary, Caigary, Alberta, Canada, Sept. 1976. 24. 25. 26. 27. 28. 29, 30. 31. 32. 33. -86- J. C. Heinrich, P. S. Huyakorn, 0. C. Zienkiewicz, and A. R. Mitchell, "An 'upwind' finite element scheme for two-dimensional convective transport equations," Internat. J. Numerical Meth. in Eng. 11, 131-143 (1977). T.d.R. Hughes {ed.), Finite element methods for convection domin- ated flows, AMD Vol. 34, ASME, New York, Dec. 1979. T.J.R. Hughes, "A simple scheme for developing 'upwind' finite elements," Internat. J. Numerical Meth. in Eng. 12, 1359-1365 (1978). T.J.R. Hughes and J. Atkinson, "A variational basis for ‘upwind’ finite elements," IUTAM Symp. on Variational Methods in the Mechanics of Solids, Northwestern University, Evanston, I11., Sept. 1978. A. Brooks and T.J.R. Hughes, "Streamline-upwind/Petrov-Galerkin methods for.advection dominated flows," Third Internat. Conf. on Finite Element Methods in Fluid Flow, Banff, Canada, 1980. T. Belytschko, J. M. Kennedy, and D. F. Schoeberle, "Quasi- Eulerian finite element formulation for fluid-structure interac- tion," J. Pressure Vessel Technol. 102, 62-69 (1980). R. H. Galltagher, J. T. Oden, C. Taylor, and 0. C. Zienkiewicz (eds.), Finite Elements in Fluids 1: Viscous Flow and Hydrodynamics (John Wiley, London, 1975). R. H. Gallagher, J. T. Oden, C. Taylor, and 0. C. Zienkiewicz (eds.), Finite Elements in Fluids 2: Mathematical Foundations, Aerodynamics, and Lubrication (John Wiley, London, 1975). R. H. Gallagher, 0. C. Zienkiewicz, J. T. Oden, M. Morandi-Gechi, and C. Taylor (eds), Finite Elements in Fluids 3 (John Wiley, London, 1978). D. K. Gartling, "Finite element analysis of viscous incompressible flow," TICOM Report 74-8, University of Texas, Austin, December 1974. 34, 35. 36. 37. 38. 39. 40. 4}. 42. 43. AA. -87- D. K. Gartling, "Recent developments in the use of finite element methods in fluid dynamics” in Computing in Applied Mechanics, AMD Vol. 18 (ASME, New York, 1976), pp. 65-92. P. Gresho,. R. Lee, and R. Sani, "On the time-dependent solution of the incompressible Navier-Stokes equations in two and three dimensions," in Recent Advances in Numerical Methods in Fluids (Pineridge Press, Swansea, to appear). C. Taylor, K. Morgan, and C. A. Brebbia (eds), Numerical Methods in Laminar and Turbulent Flow (John Wiley, New York, 1978). D. K. Gartling, "Some comments on the paper by Heinrich, Huyakorn, Zienkiewicz, and Mitchell," Internat. J. Numerical Meth. in Eng. 12, 187-190 (1978). 0. C. Zienkiewicz, “Reply to comments by Gartling," Internat. J. Numerical Meth. in Eng. 12, 191 (1978). T.J.R. Hughes and A. Brooks, “A multi-dimensional upwind scheme with no crosswind diffusion," AMD 34 (ASME, New York, 1979), pp. 19-35, A. Brooks, "Finite element solution to Navier-Stokes equation," Ph.D. thesis, California Institute of Technology. In preparation. P. M. Gresho, R. L. Lee, and R. L. Sani, “Advection-dominated flows with emphasis on the consequence of mass lumping,” in Preprints of Second Internat. Symp. on Finite Element Methods in Flow Problems, S. Margherita, Italy, June 14-18, 1976. G. B. Whitham, Linear and Nonlinear Waves, Pure and Applied Mathe- matics (Wiley-Interscience, New York, 1974). J. Donea, P. Fasoli-Stella, S. Giuliani, J. P. Halleux, and A. V. Jones, "The computer code EURDYN-1M for transient dynamic fluid- structure interaction, EUR 6751 EN , Commission of the European Communities Joint Research Center, Ispra Establishment, Italy. Preprint, 1980. T. Belytschko and U. Schumann, "Fluid-structure interaction in light water reactor systems." Preprint. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. -88- T. Belytschko (ed.)}, Fluid-Structure Interaction in LWR Systems, Internat. Congress Center, Berlin (west), Germany, Aug. 20-21, 1979. T.J.R. Hughes, "Finite element method class notes, AE/AM 108", California Institute of Technology T.J.R. Hughes, "Current trends in finite element research," Feature article to appear in Applied Mechanics Reviews. T.J.R. Hughes, "Recent developments in computer methods for struc- tural analysis," Nucl. Eng. and Design 57, No. 2, May 1980. R. D. Krieg, "An efficient numerical method for time independent plasticity," SAND 77-0943, Sandia Laboratories, Albuquerque, New Mexico, 1977. J. 0. Hallquist,"NIKE 2D: An implicit, finite-deformation, finite- element code for analyzing the static and dynamic response of two- dimensional solids,” Lawrence-Livermore Laboratory Report UCRL-52678, University of California, Livermore, Ca., March 3, 1979. R. M. McMeeking and J. R. Rice, “Finite element formulation for problems of large elastic-plastic reformation," Internat. J. Solids and Structures 11, 601-616 (1975). B. Noble, Applied Linear Algebra (Prentice-Hall, Englewood Cliffs, N.J., 1969). H. L. Schreyer, R. T. Kulak, and J. M. Kramer, "Accurate numerical solutions for elastic-plastic models," J. Pressure Vessel Technol. 101, 226-234, August 1979. T.J.R. Hughes and J. Winget, "Finite rotation effects in numerical integration of rate constitutive equations arising in large- deformation analysis," Internat. J. Numerical Meth. in Eng., in press. H. Armen, “Assumptions, models, and computational methods for plas- ticity," Computers and Structures 10, 161-174 (1979). R. D. Krieg and D. B. Krieg, “Accuracies of numerical solution method for the elastic-perfectly plastic model," ASME J. Pressure Vessel Technol. 99, 510-515 (1977). 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. -89- S. W. Key, "Hondo, a finite element computer program for large deformation dynamic response of axisymmetric solids," Sandia Labo- ratories, Albuquerque, N.M., SLA-74-0039, April 1974. D. P. Mondkar and G. H. Powell, "Static and dynamic analysis of nonlinear structures," ANSR-I, Report EERC-75-10, University of California, Berkeley, 1975. J. K. Knowles, E. Sternberg, "Notes on mathematical theory of elas- ticity, AM135 and 136," California Institute of Technology. R. D. Krieg and S. W. Key, “Implementation of a time independent plasticity theory into structural computer programs," Constitutive Equations in Viscoplasticity: Computational and Engineering Aspects," AMD 20 (ASME, New York, 1976). T.J.R. Hughes, "Stability of one-step methods in transient nonlinear heat conduction," Trans. 4th Internat. Conf. on Structural Mechanics in Reactor Technology, San Francisco, CA, August 1977. M. A. Haroun, Dynamic analysis of liquid storage tanks, Ph.D. thesis, California Institute of Technology, December 1979. T. Belytschko and R. Mullen, "Mesh partitions of explicit-implicit time integration," in Formulations and Computational Algorithms in Finite Element Analysis, K. J. Bathe, J. T. Oden, and W. Wunderlich (eds.) (M.I.T. Press, Cambridge, Mass., 1977). T. Belytschko and R. Mullen, "Stability of explicit-implicit mesh partitions in time integration," Internat. J. Numerical Meth. in Eng. 12, 1575-1586 (1978). T.J.R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient analysis: stability theory," J. Appl. Mech. 45, 371-374 (1978). T.J.R. Hughes and W. K. Liu, "Implicit-explicit finite elements in transient analysis: implementation and numerical examples," J. Appl. Mech. 45, 375-378 (1978). C. A. Felippa, K. C. Park, and H. C. Yee, “Synthesis of staggered solution procedures for coupled field analysis," 2nd Internat. Conf. 68. 69. 70. 7). 72. 73. 74, 75. 76. -90- Appl. Numerical Modeling, Polytechnical University of Madrid, Madrid, Sept. 11-15, 1978. K. C. Park, C. A. Felippa, and J. A. DeRuntz,"Stabilization of stag- gered solution procedures for fluid-structure interaction analysis," in Computational Methods for Fluid-Structure Interaction Problems T. Belytschko and T. L. Geero (eds), Applied Mechanics Symposia Series (ASME, New York, 1977). K. C. Park, "Partitioned transient analysis procedure for coupled- field problems," presented at Second Conf. on Numerical Methods in Nonlinear Mechanics, TICOM, University of Texas, Austin, March, 1979. K. C. Park and C. A. Felippa, "Direct time integration methods in nonlinear structural dynamics," Computer Meth. in Appl. Mech. and Eng. 17/18, 277-313 (1979). K. C. Park and P. G. Underwood, "A variable-step central difference method for structural dynamics analysis - Part I: Theoretical aspects," ASME Paper No. 79-PVP-120, presented at the Pressure Vessels and Piping Conf., San Francisco, CA, June 25-29, 1979. T.J.R. Hughes, K. S. Pister, and R. L. Taylor, "Implicit-explicit finite elements in nonlinear transient analysis," Computer Meth. in Appl. Mech. and Eng. 17/18, 159-182 (1979). T.J.R. Hughes and R. A. Stephenson, "Convergence of implicit- explicit algorithms in nonlinear transient analysis," to appear in Internat. J. Eng. Science. T. Belytschko, R. L. Chiapatta, and H. D. Bartel, “Efficient large scale nonlinear-transient analysis by finite elements," Internat. J. Numerical Meth. in Eng. 10, 579-596 (1976). T. Belytschko, J. R. Osias, P. V. Marcal (eds.), Finite element analysis of transient nonlinear structural behavior, AMD 14 (ASME, New York, August 1975). T. Belytschko, "Transient analysis,” in Structural Mechanics Computer Programs (University Press of Virginia, Charlottesville, VA, 1974), 77. 78. 79. 80. 81. 82. 83. 84. -9]- pp. 255-276. H. M. Hilber, T.J.R. Hughes, "Collocation, dissipation, and ‘over- shoot! for time integration schemes in structural dynamics," Earthquake Eng. Struct. Dynamics 6, 99-118 (1978). T.J.R. Hughes, T. K. Caughey, and W. K. Liu, "Finite element meth- ods for nonlinear elastodynamics which conserve energy, J. Appl. Mech. 45, 366-370 (1978). T.J.R. Hughes, “Implicit-explicit finite element techniques for symmetric and nonsymmetric systems," Proc. Conf. on Nonlinear Problems in Mechanics, University of Swansea, Swansea, U.K., Sept. 2-5, 1980. T. Belytschko, H. J. Yen, and R. Mullen, "Mixed methods for time integration," Computer Meth. in Appl. Mech. and Eng. 17/18, 259- 275 (1979). P. G. Underwood and K. C. Park, "A variable-step central difference method for structural dynamics analysis - Part II. Implementation and performance evaluation," ASME Paper No. 79-PVP-121, presented at the Pressure Vessels and Piping Conf., San Francisco, CA, June 25-29, 1979. J. P. Wright, "Mixed time integration schemes," Computers and Structures 10, 235-238 (1979). J. Isenberg, D. K. Vaughan, I., and I. Sandier, "Nonlinear soil- structure interaction," EPRI NP-945, Project 810-2, Final Report, Dec. 1978. N. M. Newmark, "A method of computation for structural dynamics," J. Eng. Mech. Div. ASCE, 67-94 (1959). -92- Chapter 4 NONLINEAR FINITE ELEMENT ANALYSIS OF SHELLS I. Introduction Within recent years, shell structures have assumed increased im- portance due to new developments in industry; in particular in off- shore structures, nuclear reactors, cooling towers, and oi] tanks. The demand for "better" design of a shell structure requires improved shell analysis procedures, and in many cases the safety of shell structures cannot be estimated unless a nonlinear analysis is undertaken. As a result of the increased speed and storage capacity of present computers, many nonlinear numerical analyses of shells have been at- tempted. Finite element procedures for plate and shell analysis now span over 20 years. A history of finite element shell analysis is traced in [1,2]. Descriptions and references to many important works may be obtained by consulting standard texts [3,4], the review articles [5], and thesis [6]. The innumerable proposed shell and plate elements demonstrate the wide dissatisfaction with available methodology. There are two classical approaches in finite element analysis of shells: (1) the direct approach whichis deduced from a classical shell theory [7]; (2) the degeneration approach [4,8] which is deduced from the field equations of the three-dimensional theory with various assump- tions. The associated shell differential equations derived from the direct approach, often higher than second order, may require a higher- order interpolation function by the finite element process. Another disadvantage of this approach is the unavailability of a convenient, -93- general enough nonlinear shel] theory. It is concluded that to develop a general enough nonlinear shell theory which can accommodate finite strains and finite rotations, and permits use of arbitrary, three- dimensional nonlinear constitutive equations, one must adopt the "degeneration" approach. Mindlin theory [9] is the most widely adopted basis of the degeneration approach, since generalization to the fully- general nonlinear case is straightforward, and an analogous two-dimen- sional formulation which includes several special cases of practical interest, namely the axisymmetric case for shells of revolution, the plane strain case for long tubes, and the plane stress case for rings and two-dimensional beam/frame structures, can be deduced fromthe three- dimensional formulation. This theory was originally employed in the linear case by Ahmad et al. [10], later by Hughes and coworkers [11-14]. Recent works in the nonlinear area which may be mentionedare [8,15-20, 39]. However, there are still certain numerical difficulties remaining in the development of a fully-general nonlinear shell theory. As has been mentioned by Hughes and Liu [20], they are (1) selective integration techniques, (2) numerical integration of constitutive equations, and (3) rank-deficiency problems associated with the in-plane rotational mode. Of course, there are other approaches to develop nonlinear shell theory. + The interested reader may consult [6,21-23] and references therein. The equations of motion derived in Chapter 2 can be used directly in shell analysis. The Lagrangian formulation will be adopted (j.e., R, = R,). For discussion purposes we will denote x as an undeformed -94- position vector and y as the deformed position vector. In Section II, the geometrical and kinematical behavior of a typical shell element will be described. The construction of "lamina" and "fiber" coordinate systems is derived along with the corresponding transformation arrays which relate these systems to the global frame. In Section III, numerical integration of the constitutive equations and the "corotational approximation" [22,24,25] are discussed. The definitions of element arrays, selective/reduced integration tech- niques, numerical integration of the lumped mass matrix and shear correction factors are discussed in Section IV. In Section V, we discuss the reduced integration Lagrange elements and the heterosis element. Sample problems involving a number of the elements are con- tained in Section VI. Il. Geometric and Kinematic Descriptions II-A. Geometry We denote the position vector of a generic point of the shell in the undeformed configuration, which is also the reference configuration, by x. It is equal to the sum of the position vector of a point in the reference surface x and a position vector based at a point in the reference surface which defines the "fiber direction" through the point X. Following the finite element discretization procedures described in Section IV-A of Chapter 3, a smooth mapping of the biunit cube into the physical shel] domain (see Figure IIA-1) is defined by the follow- I-v big UIDWOP juawala {jays |DdISAUd 3 aqno jiuniq -~96- ing equations: x(E.n.c) = x(Esn) + X(Esn.z) (1IA.1) x(En) = ay N,(En) x, (IIA. 2) Men X(Esn.t) = ) NA(Esn) X4(2) (IIA.3) x J x X,(2) = z,(c) K. (no sum) (IIA.4) z(t) = Ny(z) 2, + N(z) 2, (IIA.5) N(c) = # (1+), N(¢) = 3 (1-0) (I1A.6) In equations (IIA.2) - (IIA-6), Xo is the position vector of nodal point a; No denotes a two-dimensional shape function associated with node a; Non is the number of element nodes; x, is a unit vector emanating from node a in the fiber direction; and Z, is a "thickness function,” asso- ciated with node a, which is defined by the location of the reference surface. A lamina surface is defined by equation (IIA.1) for ¢ fixed; and a fiber line is defined by the same equation for &,n fixed. In general, the fibers are not perpendicular to the laminae. It has been found convenient to take as input the coordinates of the top and bottom surfaces of the shell along each nodal fiber (x) and x respectively) and a parameter TE[-1,+1] which defines the location of the reference surface. For example, if C=-1,0,+1 (respectively), then the refenence surface is taken to be the bottom, middle, top (respectively) of the shell. From these input data, we may calculate Xa» X Za and Z5 (a=1,2,°++sMey) from: -97- rs rn rr xX, ~ 2 (1-¢) x, + (1+ 2%) Xa (IIA.7) a _ + ~ + - Kos (xg> Xa) / Mixa > Xai (1TA.8) +_ 1 a z, = (1-2) IIx, - x, (IIA.9) - 1 +. z= - a (142) IIx, - x5 (IIA.10) . . 1/2 -|[ denotes the Euclidean norm (j.e., {I x|| = (x2 + x5 C+ x5) ). An illustration of these ideas is presented in Figure IIA-2. By virtue of the fact that the top- and bottom-surface coordinates are uniquely defined at element interfaces, there are no gaps or overlaps along ele- ment boundaries. Similarly, the position vector y of the same generic point of the shell in the deformed configuration is defined by: y(Esn.t) = y(Esn) + ¥(Esnsz) (IIA.11) n y(n) = 2 N(Esn) Ya (IIA.12) Men Y(Gsnae) =] N(Eon) Yale) (114.13) a= Y, 63) = z,(%) Y. (no sum) (IIA.14) where z,(%) and N,,N_ are defined by (IIA.5) and (1IA.6), respectively. We observe from equations (IIA.4) and (IIA.14) that the nodal fibers are assumed to be inextensible, i.e., they may rotate, but cannot stretch or contract. Therefore Y. is a unit vector. ~98- (4°3) X= (7U'F) KX 2-V D Ol4 aoDJINs WoO}oO — BIDJuNS aduaIajaI SS (92 (1+ "L'3) xX adDjuNs do} fan -99- II-B. Kinematics We define the displacement u of a generic point to be the dif- ference between the deformed and the reference configurations of the same point. It can be written in terms of the sum of the displacement of that point on the reference surface u, and the"fiber displacement" U. Therefore (see Figure IIB-1), ~ ury-xeutl (11B.1) u =y-x (1IB.2) = % - % (118.3) U =Y- xX | (11B.4) U,=Y,- X, 72 Qu, (no sum) (11B.5) b= Ya> hy (118.6) By virtue of the definition of Us the kinematics of the shell element are defined by invoking the isoparametric hypothesis that the same ex- pressions are used for kinematics as for geometry with displacement variables in place of coordinate variables. In order to define y, a trial value of U. is calculated (i.e., gtrtal) and projected radially to maintain inextensibility. The steps in the procedure are as follows (see Figure IIB-2): 9 $ atrial ° “trial Ya (X + Us y/ IX, * Us | (1B.7) U, =Y- X (IB.8) We have found the present scheme to be both economical and effective UOI}D4NBIJUod > PSWAOJIP ul QODJINS sdUaIaJaU 100- \ \ _”" I-@ D614 n een ee uolpounBijuoo pawOosepunN ul JIDJUNS adUaJaJaI (UOI}D|SUD4} |a|JOIDd) N I>< -101- unit sphere about node a Fig. I B-2 Nodal fiber inextensibility condition maintained by radial return normalization -102- in practice. II-C. Lamina Coordinate Systems In order to invoke the plane stress hypothesis, at each integration point in the element a coordinate system is defined as follows: RL ey =y e/Ily ell (IIC.1) R R R e4 = epxy ,/I epxy al (1IC.2) R_ 2k L ep = eg X ey (11C.3) where a comma is used to denote partial differentiation (e.g., Yor = ay/9é) and "x" denotes the cross product. We call this coordinate system the lamina coordinate system (see Fig. IIC-1). The e; direction is used for the purpose of invoking the zero normal stress hypothesis. We should note that e3 in general is not parallel to the ¢ axis (see Figure IIC-2), that is, the fiber direction. In order to have a global finite element formulation, an orthogonal transformation matrix q is ~ defined: q = [ey es es]! global + lamina (I1C.4) II-D. Fiber Coordinate Systems Due to the plane stress hypothesis, the rank of the material re- sponse matrix C is five, even though it is a 6x6 matrix in the gen- eral three-dimensional constitutive theory. In order to remove -103- ‘|SUOD Wa{Sks a}DUIP4009 DulluD] |DDIdA| |-5 TT ‘bi4 ¢ 7° SuOD -& (QODjJANS *{SUOD = 9) DUIWD| |DDIdA\ -104- Jaq} D Buojo swajshs ayDulpsood pulwo7 2-9 TT “bl4 aodDjINs *jsuod - U se ys: C. (A) 29 (uoysauip saqy) Sy WOAH 3 T -105- the in-plane twist mode, we define a unique local Cartesian coordinate system at each node which is used as a reference frame for rotation in- crements. The algorithm is as follows: 1. Let a, = iY. | , 1=1,2,3 (IID.1) 2. j=] (IID.2) 3. If a, > a3, then a, = a); and j = 2. (11D.3) 4. If a, > a3; j = 3. (IID.4) f 9 5. e, = ¥ (11D.5) f _ 7% n 6. e, = (Yxe,)/||¥xe;| (11D.6) f_ofio¢ 7. e, = @,xY (11D.7) This orthonormal fiber basis is chosen in such a way that if Y is "close" to the third axis of the global system, then er ep. ef will be "close" to the global axes (see Figure IID-1). With these definitions, an orthogonal transformation matrix r is defined such that for each time step the three incremental fiber displacements in the global coordinate system can be contracted to two incremental rotations. Let AU, AU,» Au, be the global incremental vector components of Ab, also let As, and AG, be the rotation increments about the basis vectors ef and er» respectively. The sign convention is defined in Figure IID-]. The relationship between AU and A@ is: (11D.8) > 2e> iH] | > @@m where r: fiber > global (IID.9) ~ -106- [fiber / e A ° e,-Y » Co A@ e f Aé —_ / €2 OF} e, / / Fig. I D-! Nodal fiber basis and -] 0 (IID.10) hom | nH Oo ' wore AD AU (IID.11) A8 = (IID.12) Equations (IID.8) through (IID.12) enable the reduction of the nodal degrees-of-freedom from six to five in the matrix incremental equilib- rium equations. This obviates the need to develop artificial in-plane torsional stiffnesses to numerically stabilize rotation about the fiber direction [17]. When beam-type stiffeners are assembled with shell elements, the bending stiffness of the beam naturally provides the in-plane torsional stiffness of the composi te structure, equations (IID.8) through (I1D.12) can be ignored, and six degrees-of-freedom are retained in the global incremental equilibrium equations. III. Integration of Constitutive Equations III-A. Numerical Algorithms The numerical algorithms described in Section III of the last chapter cannot be used directly, due to the stress components being referred to a rotating basis, Also, the zero normal-stress constraint needs to be enforced with respect to the third direction of the lamina -108- basis. So the integration of the rate type constitutive equations re- quires special considerations. We define y and w according to equations (IIA.1) through (IIIA. of Chapter 3. The orientations of the Jamina bases at time th and t, are defined, respectively by: q™: global + lamina at tr (IIIA. ~ n+1, . q : global + lamina at t,, (IIIA. ~ Consequently, the incremental transformation between lamina bases is given by: _— ntl, nyt, . . Aq =q (q')': lamina at t, > lamina at t.) (IIIA. ~ With these definitions, the numerical algorithms are: aT og ot Rt (IIIA. ctl. Mls ay (IIIA, where R = (1+(1-du)7'w) dg (IIIA. btis = Cask Yee (IIIA. and ont] = = ¥33 = (33 ~ ) Ssaig%ag) / 3333 (IITA. ij#33 The matrix R is in general an approximation to the rotation of 6) +] 1) 2) the material particle relative to the lamina basis; (IIIA.8) is used to insure satisfaction of the zero normal-stress condition in the lamina ntl = Q), basis at tae] (i.e., 733 In many situations of practical in- terest, the particle rotation will be closed to the rotation of the -109- lamina basis and thus R fs I. This is the case for what is classically known as the "small strain-finite rotation approximation. In this case, the calculation of R may be omitted. This philosophy has been exploited in the works of Argyris [21,23,26,27] and Belytschko [22,24, 25]. IV. Element Arrays IV-A. Material Tangent Matrix The element tangent matrix D is defined by equations (IVB.13) of Chapter 2. For application to shell analysis, the D matrix needs to be modified to account for the zero normal-stress condition. We define a zero normal-stress projection operator P such that the row and column of D corresponding to ¥33 are removed, 1.€., b= ppp (IVA.1) where Po p=| ~ ~63 (1VA.2) O56 13 1 0 0 oO 0) 0 1 0 oO 0 0 o 1 0 oO P = (IVA. 3) Py Pp Pz Pe Pg 0 0 0 71 0 po 0 60 OT and (IVA.4) Pr =~ ©3355 / &3333 -110- IV-B.. Strain-Displacement Matrix In order to prevent "mesh-locking" phenomena [11-14, 17-20], special treatment needs to be given to transverse shear terms. We will employ the same idea, namely reduced/selective integration techniques, which have been discussed in Chapter 3. We implement the selective integration procedure by a simple modification of the strain-displace- ment matrix. The definition of the strain-displacement matrix adopted is given as follows: B= [B).Bys--+.8, J (1VB.1) Ba 7a = ~ ba a= 1,2,+++ Ng, (IVB.2) a) Ba (~ = B, 0 0 +B, Oo 0 l Bp 8, 9 | Bs B, 0 0 B 0 ;0 B oO B= noe anne a an (1VB.3) _ — | — — 0B, Bo Be Bs (B, 0 B, 1B 0 By J - B, -B, O ' By -B, 0 ! Bo =| 0 Bz -By 10 Be -B, (1VB.4) ' -B; 0 By ,-Be 0 By | Nai ’ 1=1,2,3 B. = (1VB.5) (NaZq) (1-3) > T7498 -111- Note that Nai = aN, /2y; where y; are the components of the coordinates in lamina systems. The definition of the B,'s is found in equation (IA.8) of Chapter 3. For details consult [20,28]. For the two- dimensional case, uniform reduced-integration Lagrange elements are em- ployed. In this case,B,'s are taken to be equal to B.'s. There appears to be no general nonlinear theory which includes shear correction effects. An ad hoc procedure, which amounts to replac- ing each B. in (IVB.3) by Kl /2 B. where « is the shear correction factor, is employed herein. (For the example problems in Section VI we have em- ployed k=5/6.) The strain increments may be computed as follows: Yvec = ii Bylo (1vB.6) wvec a=l ~@ Aus where y)) f mn | Y2 er12 Yee" 4 Y3 \=( Yoo ) (1VB.7) Y5 2753 U6 J 2731 J IV-C. Transformation Matrices Since we employ a global finite element formulation, all the ar- rays have to be transformed to the global system before they can be assembled by the standard finite element procedures. The stiffness matrices and internal force vectors are defined with respect to the lamina coordinate systems. Two transformation matrices are required: -112- (— q 035 6 dof | 933. 9 Q. = { (IVC.1) ~ ~3 5 dof 5 Z Lb ~33 ~a -1 0 roer,| 0-1 (IVC. 2) 0 - I¢ 6 dof Ss =o (1vC.3) Sa = Tz 035 7 5 dof a 025 S5 ret S.=S,{| 0 -1 (Ivc.4) Lo oO IV-D. Tangent Stiffness Matrix and Internal Force Vector The element stiffness and internal force are defined as follows: k= Dhl 8 HE = 1 (190.1) +] k= TT 2 Keb | fg Q. BS D BQ, J dz dO (IVD. 2) oO -! + fint . - | | ae "eed dz dO (IVD.3) oO -l where -113- +] +] | ees dO= | | eee d& dy (laminar integral) (IVD. 4) 5 -1 -1 ~ 7 Me Yen N,e | Y3,€ ¥3,n ¥3,0 | "1 72 "12 Tvec ( 13 b= ( "22 ) a) 723 er, 31) In most situations the variation of Q. with 5 will be insignifi- cant. For this reason we can take Q. outside the fiber integral, viz, ~ c+] Oo L-l int _ Y\T . fa ~ | Qs | (B,) Tvec J d,) dO (1VD.8) 0D _-l IV-E. External Force Vector We allow for both body and surface force vectors. Body force The element body force vector is given by: body - {00%} (1VE.1) ~ -114- +] body _ oT T . ta 7 Sa | | Na b Podg de do (IVE.2) oO -! where Po is the mass density in the initial configuration, b is the pre- scribed body force vector (per unit mass), VE “Ian “Tz jy = | Xone Xan *2ar (IVE.3) 3,6 “3 “3,2 and N, 0 0 ! Nz, 0 0 NJ =] 0 Ng 0 4 0 Ny ? (IVE.4) 0 oO N, + 0 O Naz, Surface force The element surface vector is defined by fut pgsurhy (IVE.5) +] top surf = si | Nihj, dd, c= (1VE.6) ~ ~ “ -1 bottom oO j, = lly-xy yl surface Jacobian (IVE.7) where h is the surface force vector per unit surface area. Since we consider rather general nonlinear behavior, the fiber integrals need to be evaluated by a numerical integra- tion technique instead of analytically. When the shell consists of one homogeneous layer, Gaussian quadrature is most efficient. If it is desired to include the outermost fiber points (i.e., c=+1) then -115- Lobatto rules are most accurate. If the shell is built up from a series of layers of different materials, then Gaussian rules may be effectively used over each layer. If there are a large number of approximately equal-sized layers, then midpoint rule on each layer should suffice. IV-F. Stress Resultants Bending moments, membrane forces, and transverse shear resultants may be computed at the lamina quadrature points of the "normal" rule at which stresses are stored. Let Ee (Ey Ep) be the quadrature points, 1< 2 the resultants are: <nea: =~ = “int? moments +] myglEq) = | tyglEget) 2(Egst)de 2 (Eg)s TsasB<2 (IVF.1) -| : 2(Ep st) = Ny(t) 2°(E,) + N(c) z (Ey) (IVF.2) te ye Mey ot z (Eo) = es N,(E,) z, (IVF.3) z (Ep) = z"(Eg) - 2 (E,)1/2 (IVF.4) membrance forces +] nyglEp) = | tyg(Eget) dc 2, -(Ep)» 1sa,e<2 (IVES) shears +] a t) dt z (E,) (IVF.6) ack as‘ 2k’ as 2 . The sign conventions for the stress resultants are illustrated in -116- SJUD}|NSa4 SSB1jS JO} SUOYU@AUOD UBIS |-4 AT ‘B14 -117- Figure IVF-1. Iv-G. Mass Matrix If the consistent mass matrixis desired, Gaussian integration rules are used. The consistent mass matrix is: +] Map = | | NAN, Jy de dO (1VG.1) co - where j. and N. are defined by (IVE.3) and (IVE.4), respectively. 0 a If lumped mass is required, Lobatto rules are employed. However, the At restriction is stringent for explicit transient calculations [13,38]. Therefore, the following algorithm is used to obtain a lumped mass matrix for shells such that a much larger time step can be employed with no convergence rate loss. Algorithm rot a) Define variables; M, translational mass = rotational mass, yaise - 1. Clear necessary variables and arrays. 2. Loop on laminar integration points. 3. Loop on fiber integration points. rot _ yrot 2.77 ~ % 4, M. M, + PN, JglEgoNp op) Wo. where Wo is the weight of the gth 5. MCISP — yrot a a integration point. if not heterosis. disp _ ,disp 2.47 ~ % 6. Ma Mo + Pg (N, +P. (0)Ng) Jo lEg Mp Fp) We if heterosis. : _ disp _ 7. If heterosis and a=9, M. Po, (Here P's are the serendipity shape functions. ) b) Normalization 1. rot. yen yrot a=1 ? -118- n AG en * 2. yaisp _ yetse a=] rot rot, XM 3. Ma © Ma [erot! disp disp, XM 4. M * Ma (disp! where XM is the total mass. c) Adjustment to rotational mass: ] + _l -\2, 1 7+ +)2 la, *q (z+ 23) + 72 (Zz, - 4) - Men 4 2. Zettective — L, (2, - 24) / Men 3. A = effective area = volume / ~ effective 4.0, = max{c. 3 A/8} rot rot 5, M. + M. O., For further details of mass lumping for shells and the critical time step calculation, see [13,38]. VY. Shell Elements The elements employed in this work are generalizations of ones that have been proposed by Hughes [11-14]. IV-A. Reduced Integration Lagrange Elements and Heterosis Elements The lamina shape functions and quadrature rules for the Lagrange elements are shown in Figure VA-1. In each case the appropriate re- duced rule is one order lower than the normal rule. The corresponding lamina shape functions and quadrature rules for the two-dimensional cases are shown in Figure VA-2. For this case uniform reduced inte~ gration is optimal. As has been discussed in detail by Hughes and -119- Sjuawaja jjaus sHunsbo7 j-yv A Bis EXE oXd jn UDISSNDO paonpa DX exe oxed a|nd UDISSNDS |OWOU o1qnoiq 21,Duponbiq ADaUl}!q SUOILOUNY adpus DuUIWD| ~120- SJUBWA|a |]aYS |DUO!SUSWIP-OM] Uolj}OJBajyu! paonpas wuosiuN 2-y A bl4 ain ulod— ulod — ulod — yu! £ mI ¢ a | auntosponb suoljoun} 91anod d1LD4aponb ADaul _ IQ 1}D4p I] adpys -121- coworkers [11-13,20], selective reduced integration elements behave well in thin shell applications, but may occasionally engender rank defici- ency [11-13]. Research has been undertaken to efficiently remove the mechanisms. Hughes and Cohen [11,12] proposed the heterosis concept in which the 9-node Lagrange shape functions are used for rotation fields and the 8-node serendipity shape functions are used for the displace- ment fields. The resulting element possesses correct rank and behaves well in the thin shell limit. The heterosis element is implemented by first constructing the arrays for the selectively integrated Lagrange element (say, @.g.; Kiag? Atiag? etc.), then a projection matrix, H,is constructed from the serendipity shape functions associated with the element-boundary nodes. For details of H see [11,12]. The heterosis arrays are then defined by _ yl khet ~ 4 Klag H (VA.1) T thet ~ ATiag VI. Sample Problems All computations were performed on an IBM 3032 computer at the California Institute of Technology Computing Center in double precision (64 bits per floating-point word). If linear elastic properties were used, two-point Gaussian fiber integration was employed. If plasticity theory was used, four-point Gaussian fiber integration was employed. The following notations will be used for identifying the types of ele- ments in three dimensions: -122- Ul - 4-node, uniform reduced integration S1 - 4-node, selective reduced integration S2 - 9-node, selective reduced integration HS2 - 9-node, selective reduced integration heterosis elements The full definition of the rotation matrix R and the corotational approximation R=I were both employed. In no case were the differ- ences discernible. The material model is considered to be homogeneous and linearly elastic unless otherwise specified. The shear correction factor, «, was taken to be 5/6 throughout, and for three-dimensiona] case, the five-degree-of-freedom nodal system was employed; for the two-dimensional case the three-degree-of-freedom nodal system was employed. For dynamic analysis, the Newmark parameters were taken to be 8 = 1/4 and y = 1/2 throughout, unless otherwise indicated. VI-A. The Elastica The post-buckling behavior of a clamped beam column was analyzed with a mesh of twenty 2-node elements. The analytical solution is available in Timoshenko and Gere [29]. The beam column is tilted ini- tially off center with a slope of 1:500. This imperfection is intro- duced to initiate a non-trivial solution after bifurcation (buckling) . _ Data employed in the calculations are given as fol- lows: Young's modulus = 4.8x 10°; Poisson's ratio = 0, length (L)=10, and thickness (h) = 0.01. Load-deflection results and deformed pro- files are shown in Figure VI-A. DIIISDje9 BY, VY TA ‘Bly S8}ifOud Paw4djap uojdajjap sa poo} ‘17M VM O'| 85 9° , . Ol gs oO 2 O TS / 3 ldg7e B G2B8°1 ay | Jeo ee fe i a oO S 09S lee, ~6—-6—5--6-2- - G‘| e8c'l ae J JUaW3}a Byiuiy O'?2 a yODxa 4G°2 “J Ott j os -124- VI-B. Pure Bending of a Cantilevered Beam A thin cantilevered beam was subjected to an end moment. The analytical solution for this problem is given by the deformed midsur- face curvature = ML/EI, where M is the applied bending tip moment, L is the length of the beam, I is the moment of inertia about the axis of bending and E is Young's modulus. The end moment was increased linearly up to a final value of 27. Two-hundred load steps were em- ployed and approximately two iterations per step were needed for con- vergence. The parameters of the beam were chosen such that the exact final diameter to length ratio (D/L) equaled n! = 0.3183. Both. twenty 2-node elements and ten 3-node elements were employed to model this beam and the results were found to be virtually indistinguishable. Deformed profiles are shown in Figure VI-B. The approximate D/L equaled 0.3182. VI-C. Large Deflection of a Diamond-Shaped Frame A pinned-hinged diamond shaped frame was analyzed, and results were compared with analytical and experimental data from [30] (see Figure VI-C). Thirty 2-node elements were employed to model one- quarter of the frame structure due to the double symmetries. Seventy- six equal load steps were employed. The average number of iterations per step over the first forty steps was less than three. Subsequently the number of iterations necessary for convergence increased. During the last six steps between sixteen and thirty iterations per step were required. -125- WDA J9AB|IJUDD D JO Bulpuaq aundg SjuUdWa}a apou-¢ OC} PUD SjuaWAa|a apou-Z OZ oD iG % i a iw (i) oO (B) (7 fy) ’ ny tH) 4’. + ‘) vv Ss) gm Obl4 O'8- O'9- Ov- O'c- Oxe awWD4j Paddys-puowoIp D jo uo}JIatJap a640q 5 DR-Biy ‘1D Ja suIyuar cy -P sjuawala apou—zg WoOdj 340 SyjNSas JDINNAjDUD puo DJOP jUaUI|Jadxa Aj4siys soy sadoys pajsajjap uolsp6u0|@ jo;uOZ0Y mV ALTA vo 20 ZO 31'°0 800 9) zo T T T T fs) Iafjid =2 ul Zfci= 1 z'0 (MIOT=4) Yt avi = I "7 oO = 4 7 0°9- isd ,O1xO'E = 2 +0 . JB b o’e- Yd ® a 490 4 foe] 4 U ye wo / 380 rN] Y? 7 4 ot “ ol o’s UOIjDOS}UOD |DINjIaA WAS iA] 3¢) 270 310 80'0 (9) o'6 T T T T 0 wasudgy gv DIOP (Hjyuswjssdxao 0 Asooyy apaul] —-— . Uo}98] 9p flows ——— 420 oz! uoHoe;jep 210) —— 4 ¥'0 CY bO'sl ye , git aA oe, 490 4 «4 7 i ‘ ee +90 oe / Z O71 -127- VI-D. Clamped Spherical Cap Subjected to a Concentrated Load A thin spherical cap which is clamped around the circumferential edge is subjected to a concentrated force at the vertex. Twenty and thirty 2-node elements and ten and fifteen 3-node elements were employed for the spatial discretizations. Deformed profiles and load-deflection results are shown in Figure VI-D. The different spatial discretizations did not produce discernible difference in results. However, refining the load discretization from fifty to one-hundred did bring the results in closer agreement with results of other investigators [31-33]. VI-E. Clamped Circular Plate Subjected to Uniform Pressure A clamped circular plate subjected to uniform pressure (noncon- servative) was analyzed using thirty-two and sixty-four 2-node elements. The displacement versus load results and bending moment versus displacement results are compared with analytical results [34] in Figure VI-E. Twenty-four equal load steps were used. Both meshes gave virtu- ally identical results. VI-F. Plate Strip under Uniform Load This problem is actually one-dimensional and an analytical solu- tion is available in Timoshenko and Woinowsky-Krieger [34]. We applied the plane strain assumption to a mesh of five 4-node elements for half of the plate. The plate is simply supported. Pertinent data and results are shown in Figure VI-F. The load was treated as nonconservative. Comparison is made with the analytical solution and resuits of Horrigmoe [6] who employed hybrid elements. As expected, both S1 and R = 4.76in h = 0.01576 in H = 0.0859 in @ = 10.9° E = 10 x10® Ib/in? y=03 deformed profiles 100 T T T * Bathe — — Stricklin (thirty shell elements) f . @ Mescall (finite difference solution) 80 ©. present study 100 steps xd x present study 50 steps / (twenty and thirty 2-node xd / elements, and ten and / fifteen 3-node elements 60/- yield the same results) xf 2 / = x ou . linear 40; solution a 7 iv ¥ 20 7 Oo 1 | [ i 0.4 0.8 1.2 1.6 2.0 Wo /H Fig.YVI D Clamped spherical cap subjected to concentrated load -129- aunssaid Wiojjun 0} pajoafqns ayD|d spjnd41I9 padwd|d J Mm bly juawadDidsip sA jJuawow Bulpuagq 8'I y/ Om 91 Zl 80 vO O 4iDauljuoU ‘JuUaWOW jajuao Jpau!| | l l I O / 4Daul| ‘juawow abpa x 4Daul|uoU ‘juawow abpea PpDo| SA juawWadD;|dsIp (,4a)/A,0D pool Zi ol 8 9 + 2 0 l T T T T l O - 3¢7°0 = 4 8° Apau||uoU sjuawa|a Og apou -z nN = ANQ-AXIS = sy > pun omj-Ayalyy 4JD9UI} *‘Apnys juasaud @ | 7 Jabaisy-Aysmoulom =] ! pup oyuaysow! | — | | | | | | O'2 ae | k g0l = 3 q R (snipo4) 0°60) . > ° tL t bi t ] -130- poo| wuoylun vspun dis aojig 4 TA Bl4 [ww] Om Ov o¢ Oz 1 ; a. Ol __ 6 NS uounies 4Daul| (IN pud |S) widaq jO }/04 104 | = SJUBWBJB BPOU-f BAI} Jo Apnys juasaud [D914A{OUD (49'tyy) WW *z]=y WW O°'BOG =) aoub104—— cosa i guw/NxG8'90z¢=3 Oe - 4¢--—~ nf ~~ >~4 oe b l Ov [zuu/n] b -131- U] elements yielded virtually identical results for this case, even though the S? elements took approximately two iterations per Toad step, whereas the Ul] elements took 4 to 5. Nevertheless, CPU time for the S] elements was more than double that for the U1] elements. For small prob- lems of this type, all other things being equal, the CPU time is dominated by the number of numerical integration points. VI-G. Shallow Circular Arch under Concentrated Load The problem statement of this one-dimensional problem is depicted in Figure VI-G. This problem has been solved by Bathe et al. [35] using 8-node continuum elements. We employed sixteen 4-node elements for half of the arch. Aplane strain assumption in the out-of-plane direction was used. As in the previous analysis, the results for S1 and U1 elements were virtually identical. In each case snap-through occurred at P = 35.8 1b. The number of iterations per step was 1 to 2 for the S] elements and averaged slightly over 4 per step for the U1 elements. Again, however, the $1 analysis required more than double the CPU time of the Ul] analysis. The load-deflection curve was compared with the results due to Bathe et al. [35] and was shown in Figure VI-G. VI-H. Clamped Square Plate under Uniform Load The problem configuration is defined in Figure VI-H. Plotted results are for four HS2 elements and sixteen UI elements. We also used meshes of four $2 elements, and sixteen and sixty-four S1 ele- ments. In the latter three cases, the results fell pointwise between the former two plotted cases. For clarity, the latter three cases were omitted from the plotted results. Our finite element results were -132- POO] Pa}OsjusdU09 JapuN Yop JO|NIUID MO|;DYUS 9 TL'Bl4 [ut] Om T | ! ° ULGO'I=H Ul Ope = 7] co =4 (yidap) Ul Q’|=q isd Ol=3 (M9141) Ul gi/e =y ol oJ6EC 229 Ul piles) =u 402 4 y A UOS|IIM PUD WWIDY ‘ayjog “, ee (IN pun Is) --- ==> youd 40 {DY 40j yO baun SJUIWI]9 BPOU-p UBAa{XIS Apnis juasaud l Ov [Al] d - 133- POO; W4OsluNn sapuN ayojd aipnbs padwojy W TA ‘Bly [ww] yYaWaoD|dsIp 1a}Ua9 Ov ox O'? 01 0'0.. I I T exe y IN y N A b t CY lo = LLL¢ LL LLL zL Wi 42°0 4 1 / Ye y / “Uf y oF hes J f ip a, f pe 4 4 1 rd y : - + PIIVITT TT y ‘ YOHN|OS 4Dau}| 470 | 5 i i ADM ------- ww 2 =Y Hj jj ww oool=1 ff LE co: ff? guiw/by ,O1xZ =a //s HI: iff: / : i : (S}|]NSa1 sAogo UaaMjaq a}] SJUdWI|a [S @POu-p 4NO}-AjXIS puD Uaa}xXIS ‘sjuaWa|a ZS Apou-~G uno} ‘os}o) ajOjd yO 4a,4Onb 104 S{UIWA]a |] PPOU—fp UaayxIS —o— BJO{d 4O 4a,4OND 10; S{UBWIJS SSH SISOIDJ9Y INO] «reapers ! | 490 ~134- compared to the classical Rayleigh-Ritz solution obtained by Way [36] and the finite element results of Kawai and Yoshimura [37]. For this problem we experimented with treating the load as conservative and non- conservative. This produced no discernible differences. We also ex- perimented with omitting the qT and C matrices in the calculations. The results were almost identical, even though the number of iterations per load step was increased from an average of 3 when we included both q and C, to about 7 or 8. VI-I. Hinged Cylindrical Shelli under Concentrated Load Figure VI-I shows a circular cylindrical shell with a concen- trated central load applied on the convex side. The straight edges are hinged, whereas the curved edges are completely free. Due to sym- metry, four-element meshes of S2 and HS2 elements and sixteen 4-node elements (both S1 and U1) are used to model one-quarter of the shell. All of our results were in close agreement and could not be distinguished on the scale of the plot. Comparison is made with the results obtained by Horrigmoe [6] and Bathe and Bolourchi [15]. VI-J. Transient Elastic-Plastic Response of a Simply-Supported Plate This problem was solved previously by Liu and Lin [40] who used small deflection theory, ignoring membrane effects. In our calcula- tions we assumed full geometric as well as material nonlinearity. One would anticipate some stiffening in this case due to the tensile stresses developed as the plate deforms. Since the center displacement is of the order of half the plate thickness, one would expect small, but not insignificant differences. As may be seen from the results of -135- 2540 mm 254mm l2.7 mm 0.1 rad 3102.75 N/mm2 0.3 Yr mos T- DW t tt straight edges are hinged and immovable, curved edges are free 2.5 rr ee TT rn es a Bathe and Bolourchi Pmax = 2:24 kN os >ok ON wet Z Horrigmoe 1.5e present study _ =z x oO 1.0 sixteen 4-node elements — (St and Ul) also four 9-node elements (S2) and four heterosis a 0.5 elements (HS2) Pmax 72-20 kN in each case 0.0 l | I l i i | l i 0.0 5.0 10.0 Center displacement [m m| Fig. WI | Hinged cylindrical shell under concentrated load “aydjd paysoddns—Ajdwis 0 jo asuodsas jualsunsy "PIN Bid (99S, O]) Swit -136- 8 9 b 2 0 $}UBWI|a |S UaE4xXIS I T 0 HoNdxe@ 988 9.01xGZS°S] _ \ Holdwy 988 goOlx ogz’2s . sspW padun| ‘41o1)dxa tp’ d SSDU JUSPSISUOd ‘4191;dui| —.— b ° Oe ; is oS Y= b SspwW paduwny ‘y1o1dwy ——— = pul/22®S-4) ,O1xeec’2 = 0 GLE1 UWI] PUO NIA] oveeeeeee + 3° 2 (ssa4js pyaih) isd Ol xO = ko . cox A LAY f =. Isd Ol = a : \ ye 4sD1a Jo. Oo (M214) Ugo = Y \. A DIFD4S Q uiol = 7 oiaspja NN 4 es .y 4 491 ‘ \ hj e = ‘, * - if: b I YX. me doz = INS f 5 J y SY “a 4o°2 y, g oo 2 d4so}d 4 eae see Alyoajpaad—di4spja y +48°2 4) y 4 TI ISD IIA I 7 -137- Figure VI-J, the expected behavior is produced by our calculations which tend to be slightly stiffer than Liu and Lin's with respect to peak response. This is more pronounced in the elastic perfectly plastic case than in the purely elastic case as may be seen. This is also reasonable, since larger membrane effects would tend to develop due to the smaller bending resistance in plasticity. Comparisons between lumped and consistent mass matrices, and im- plicit and explicit algorithms tend to be quite good. Despite the fact that the explicit time step was taken to be one-fourth the implicit, an economy of approximately a factor of seven was noted for the explicit calculations in the elastic-plastic case. In the purely elastic case this was reduced to a factor of four favoring the explicit technique. VI-K. Impulsively Loaded Elastic-Plastic Strip Experimental results for this problem were given by Balmer and Witmer [41] and results of a finite element calculation were presented by Belytschko and Marchertas [42] who used an elastic, perfectly plastic model with artificial viscosity. In our calculations, we ex- perimented with different types of hardening, but did not include any viscous effects. Another difference between the Belytschko-Marchertas calculation and ours was transverse shear effects which were not in- cluded in theirs, but were included in both the kinematics and consti- tution of ours. A comparison of results is presented in Figure VI-K. The displacement is quite large, being approximately six times the strip thickness. Our peak responses tend to be in good agreement with the experimental value, although they tend to occur at somewhat earlier times. The closest agreement we are able to get is attained with 0.8 0.7 (in) Q.5 center displacement 0,1 0.0 f o& © Balmer and Witmer 1964 (exp’t.) [/ —-—Belytschko and Marchertos 1974 / — —elastic-perfectly plastic -“isotropic hardening present kinematic hardening study ------ combined isotropic- kinematic hardening ay “ 1 WAP ZZZZZZ ZR NX \ 0.125" = a ) r "| |p =2.61xI04Ib-sec/int Vy = 5290 in/sec 0, =4.14x109 psi At =10°6 E = 1.0410" psi ten Ul elements is y =0.3 | l | Ep=.667 X10 psi . explicit, lumped mass 0.0 0.4 0.8 1.2 time (msec) Fig.MI-K. Impulsively loaded elastic-piastic strip. I.6 -139- purely kinematic hardening, as may be seen. Considering the ambiguities in material modeling, the calculation is felt to be quite reasonable. ~140- References 1. 10. lt. S. F. Pawsey, "The analysis of moderately thick to thin shells by the finite element method," SESM Report No. 70-12, Dept. of Civil Engineering, U.C. Berkeley, August 1970. H. H. Dovey, "Extension of three-dimensional analysis to shell structures using the finite-element idealization," SESM Report No. 74-2, Dept. of Civil Engineering, U.C. Berkeley, January 1974. R. H. Gallagher, Finite-Element Analysis Fundamentals (Prentice- Hall, 1975). 0. C. Zienkiewicz, The Finite-Element Method, Third edition (McGraw-Hi11, London, 1977). R. H. Gallagher, "Shell elements," Proc. First World Congress on Finite-Element Methods in Structural Mechanics, Bournemouth, England, 1975. G. Horrigmoe, "Finite-element instability analysis of free-form shells," Report 77-2, Div. Structural Mechanics, Norwegian Insti- tute of Technology, University of Trondheim, Norway, May 1977. P. M. Naghdi, "The theory of shells and plates," Report AM-71-1/2 Div. Applied Mechanics, U.C. Berkeley, 1971. u E. Ramm, "A plate/shell element for large deflections and rotations, in Formulation and Computational Algorithms in Finite Element Analy- sis (K. J. Bathe, J. T. Oden, W. Wunderlich, eds.) (M.I.T. Press, 1977). R. D. Mindlin, "Influences of rotatory inertia and shear on flexu- ral motions of isotropic, elastic plates," J. Appl. Mech. 18, 31-38 (1951). S. Ahmed, B. M. Irons, and 0. C. Zienkiewicz, “Analysis of thick and thin shell structures by curved finite elements," Internat. J. Numerical Meth. in Eng. 2, 419-451 (1970). T.J.R. Hughes and M. Cohen, "The 'heterosis' finite element for plate bending," Computers and Structures 9, 445-450 (1978). 12. 13. 14, 15. 16. V7. 18. 19. 20. 2l. 22. -141- T.J.R. Hughes and M. Cohen, "The ‘heterosis' family of plate finite elements,” Proc. ASCE Electronic Computations Conf., St. Louis, Mo., August 6-8, 1979. T.J.R. Hughes, M. Cohen, and M. Haroun, "Reduced and selective integration techniques in the finite element analysis of plates," Nucl. Eng. and Design 46, 203-222 (1978). T.J.R. Hughes, R. L. Taylor, and W. KanoknukulChai, "A simple and efficient element for plate bending," Internat. J. Numerical Meth. in Eng. 11, 1529-1543 (1977). K, J. Bathe and S. Bolourchi, "A geometric and material nonlinear plate and shel] element," J. Computers and Structures, to appear. B. Brendel and E. Ramm, “Linear and nonlinear stability analysis of cylindrical shells," Internat. Conf. on Eng. Applications of the Finite Element Method, Hévik, Norway, May 9-11, 1979. W. Kanoknukulchai, "A large deformation formulation for shell analysis by the finite element method," Ph.D. Thesis, University of California, Berkeley, November 1978. H. Parisch, “Geometrical nonlinear analysis of shells," Computer Meth. in Appl. Mech. and Eng. 14, 159-178 (1978). H. Parisch, "A critical survey of the 9-node degenerate shell ele- ment with special emphasis on thin shell] applications and reduced integration," Computer Meth. in App}. Mech. and Eng. 20, 323-350 (1979). T.d.R. Hughes and W. K. Liu, "Nonlinear finite element analysis of shells. Part I - Three-dimensional shells; Part II - Two-dimensional shells," to appear in Computer Meth. in Appl. Mech. and Eng. J. H. Argyris, H. Balmer, J. St. Doltsinis, P. C. Dunne, M. Haase, M. Kleiber, G. A. Malejannakis, H-P. Mlejnek, M. Muller, and D. UW. Scharpf, "Finite-element method--the natural approach," Computer Meth. in Appl. Mech. and Eng. 17/18, 1-106 (1979). T. Belytschko, L. Schwer, and M. J. Klein, “Large displacement transient analysis of space frames," Internat. J. Numerical Meth. 23. 24. 25. 26. 27. 28, 29. 30. 31. 32. -142- in Eng. 11, 64-84 (1977). J. H. Argyris, P. C. Dunne, G. A. Malejannekis, and D. W. Scharpf, "On large displacement--smal] strain analysis of structures with rotational degrees of freedom," Computer Meth. in Appl. Mech. and Eng. 14, 401-451 (1978). T. Belytschko and L. Glaum, “Application of higher-order corota- tional formulation for nonlinear finite element analysis," Compu- ters and Structures 10, 175-182 (1979). T. Belytschko and B. J. Hsieh, "Nonlinear transient finite element analysis with convected coordinates," Internat. J. Numerical Meth. in Eng. 7, 255-271 (1973). J. H. Argyris, "Recent advances in matrix methods of structural analysis," Progress in Aeronautical Sciences, Vol. 4 (Pergamon Press, Oxford, 1964). J. H. Argyris, "Continua and discontinua," Proc. Conf. on Matrix Methods in Structural Mechanics, Wright-Patterson A. F. Base, Ohio, October 1965. T.J.R. Hughes, "Generalization of selective integration procedures to anisotropic and nonlinear media," Internat. J. Numerical Meth. in Eng., to appear. S. P. Timoshenko and J. M. Gere, Theory of Elastic Stability (McGraw-Hill, New York, 1961). J. A. Jenkins, T. B. Seitz, and J. S. Przemieniecki, "Large deflections of diamond-shaped frames," Internat. J. Solids and Structures 2, 591-603 (1966). K, J. Bathe, E. L. Wilson, and R. H. Iding, NONSAP--A structural analysis program for static and dynamic response of nonlinear systems, UC SESM 74-3, Dept. of Civil Eng., University of California, Berkeley, 1974. J. F. Mescall, "Large deflection of spherical shells under concen- trated loads," J. Appl. Mech. 32, 936-938 (1965). 33. 34, 35. 36. 37. 38. 39. 40. 4]. 42. -143- J. A. Stricklin, "Geometrically nonlinear static and dynamic anal- ysis of shells of revolution," High Speed Computing of Elastic Structures, Proc. IUTAM Symp., University of Liege, 1970. S. P. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd edition (McGraw-Hill, New York, 1959). K. J. Bathe, £. Ramm, and E. L. Wilson,"Finite element formulations for large displacement and large strain analysis," Report UC-SESM- 73-14, Dept. Civil Eng., University of California, Berkeley, Sept. 1973. S. Way, "Uniformally loaded, clamped, rectangular plates with large deformation," Proc. 5th Internat. Congress Appl. Mech., Cambridge, Mass., 1938. T. Kawai and N. Yoshimura, “Analysis of large deflection of plates by the finite-element method," Internat. J. Numerical Meth. in Eng. 1, 123-133 (1969). R. D. Krieg and S. W. Key, "Transient shell response by numerical time integration," Internat. J. Numerical Meth. in Eng. 7, 273-286 (1973). E.D.L. Pugh, E. Hinton, and 0. C. Zienkiewicz, "A study of quadri- lateral plate bending elements with 'reduced' integration," Inter- nat. J. Numerical Meth. in Eng. 12, 1059-1079 (1978). S. C. Liu and T. H. Lin, "Elastic-Plastic dynamic analysis of structures using known elastic solutions," Earthquake Eng. and Structural Dynamics 7, 147-159 (1979). H. A. Balmer and E. A. Witmer, “Theoretical-experimental correla- tion of large dynamic and permanent deformations of impulsively loaded simple structures,” FDL-TDR-64-108, A. F. Flight Dynamics Laboratory, Wright-Patterson A. F. Base, Ohio, July 1964. T. Belytschko and A. H. Marchertas, "Nonlinear finite element method for plates and its application to dynamic response of re- actor fuel subassemblies," J. Pressure Vessel Technol. 96, 251- 257 (1974). -144- Chapter 5 NUMERICAL EXAMPLES OF FINITE ELEMENT ANALYSIS OF INCOMPRESSIBLE VISCOUS FLOWS BY THE PENALTY FUNCTION FORMULATION I. Introduction In this chapter, sample numerical examples are presented and com- pared with available data to illustrate the accuracy and versatility of this technique, namely, finite element analysis of incompressible vis- cous flows by the penalty function formulation. As has been mentioned in Chapters 2 and 3, the selection of the penalty parameter % 7s not trivial. Clearly, } must be large enough so that the compressibility and pressure errors are negligible, yet not so large that numerical i11- conditioning ensues. Dimensional analysis reveals that, for Stokes flow, } should be picked according to the relation XA = Cy (1.1) where c is a constant which is computer word length dependent only. Numerical studies reveal that for floating-point word jengths of 60 to 64 bits an appropriate choice of c is 107. However, in the presence of the convective momentum which generally dominates the viscous term, the criterion for the Navier-Stokes equation is empirically defined as: dA = c max{u,uRe} (1.2) where yp is the dynamic viscosity and Re is the Reynolds number as given by Re = UL/v (1.3) -145- Here U and L are "characteristic" velocity and length, respectively, v = p/o is the "kinematic viscosity," and 9 is the density. Usually U was taken to be the maximum expected velocity and L the maximum dimen- sion. The solution strategy used for the steady Navier-Stokes equations is an incremental Newton-Raphson scheme in which density is used as a "load parameter," see [1-4] for details. The solution of the Stokes problem is the initial guess. Within each load level, iterations are performed until convergence is achieved, and the converged solution is used as the initial guess for the next load level. It should be noted that the matrix to be inverted [equation (IVA.12) of Chapter 3] is non- symmetric, but possesses a symmetric band-profile structure [4]. Even though it is well known that, under appropriate hypothesis, the Newton- Raphson scheme exhibits second-order convergence of the iterates, the reforms and factorizations of the nonsymmetric stiffness are very expen- sive. Furthermore, the storage requirement is extremely demanding. Research in developing more effective solution algorithms has been a very active area of late (see, e.g., [5-7]). The solution strategy we used for the transient Navier-Stokes equa- tions in this chapter is a one-step linear implicit/nonlinear explicit operator splitting algorithm. The advantages of this linear-implicit predictor-correctcr method are that the matrix is symmetric and need be formed and factorized only once; also the critical time step is indepen- dent of the Reynolds number. The computer program developed may be run at a constant (input) time step, or at a step redefined adaptively, for each ty according to the formula below: ~146- Nsd fu, | at <1 /( y =) (1.4) j=1 Ei where Ue is the velocity component in the B direction and hey is the element length in the E; direction. It is important to cut down on reform and factorization costs when At is being selected adaptively. We have employed a scheme such that, instead of refactorization, y is redefined to compensate for the step-size change. Specifically, we proceed as follows. Let Ate act and Yeact denote the values of At and y respectively, used during the last factorization, and let A (1.5) c= Yfact Ceact Based on the velocity field at step n, Atonitis calculated according to (1.4). Define = c/At (1.6) Yn+1 crit For stability, y must be between 1/2 and 1. If Ynt1 © [pl], do not n+] refactorize, but set y = Vn] and At = At in AF. If, on the other crit hand, Yq] & Epelds set At = Atinite y = 3/4 and refactorize. This value of y is picked to reduce the likelihood of refactorization in subsequent steps. There are other procedures along these lines under investigation [8-10]. There is a drawback to the present schemes due to the convection stability condition (1.4). For example, if we are interested in an essentially steady flow, and the length of time interval, T, required ~147- to attain steady conditions engenders many steps (e.g., see cavity flows in the next section), a "fully implicit," unconditionally stable scheme is no doubt superior. Gresho et al. [11] and Smith and Brebbia [12] have described such fully implicit schemes. However, we believe that the storage (due to nonsymmetric matrix) and computational effort engendered by fully implicit algorithms are prohibitive in most cases, as often accuracy dictates taking as small a time step as the convective stabil- ity condition. To solve large, implicit, three-dimensional problems, iterative or at least partially iterative methods seem to be a necessity. Although there is much interest in this topic, a generally reliable iterative scheme, competitive with direct schemes, does not yet seem to be avail- able. Attempts have been made to overcome the mentioned difficulties; see [5-7,13-18] and references therein for details. A least-squares type smoothing procedure [19] is used to perform the necessary filtering of the "checkerboard mode" [4] of the pressure field. A comprehensive study of such techniques has been performed by Lee et al. [20]. We employed the same ideas in [20] with slight modifi- cation for our purposes. As these have been documented in [4], the smoothing procedures are omitted here. II. Sample Problems All computations were performed on an IBM 3032 computer at the California Institute of Technology Computing Center in double precision (64 bits per floating-point word). The penalty parameter is picked according to (1.1) or (1.2). Variable time stepping and "“upwinding” are employed throughout the transient analysis problems, according to -148- equations (I.4) to (1.6), unless otherwise specified. Lumped mass matrices are used for two-dimensional problems, and row sum mass matrices are used for axisymmetric problems. II-A. Driven Cavity Flow Several investigators have studied the problem (see [3,4,20-23]). A problem description is shown in Figure IIA-1. We observe that the boundary conditions are discontinuous at the upper corners. Two meshes of 4-node elements, employing different approximations of the boundary conditions, are shown in Figure IIA-2. The midplane velocity profiles employing the 20x20 and 20x21 meshes for Re = 100 and 400 are presented in Figure IIA-3. As can be clearly seen from Figure ITA-3, different treatments of the boundary conditions can result in signifi- cant quantitative differences, especially as the Reynolds number is increased. By virtue of the sensitivity of the result to the treatment of the boundary conditions, care must be taken in interpretation. For this problem, the fully ‘implicit Newton-Raphson method without "upwind" technique is employed. The boundary condition problem may also create discrepancies when different elements are compared. | Consider a 10x10 mesh of square, 9-node elements. This mesh would have the same number of degrees-of- freedom as the 20x 20 mesh of 4-node elements, and it may seem appro- priate to compare results. However, one should keep in mind that the setting of the nodal boundary conditions in identical fashion actually implies different representations along the vertical edges of the corner elements (c.f., Figures IIA-2 and IIA-4) since a linear variation in velocity is employed for the 4-node element, whereas a quadratic -149- ‘uoljdizosap WlalTqorg smo Aylaeo uaAlIq °*| VIt 6 I -V—DlA me! | > O= en=!| n O=¢n O=en Q=!n O In A/\ =94 O=en‘{=!n *suoTJIpUuOD ATepunog Tauiod snonuiucsstp jo suotjeurtxoidde jyuarazzIp sutsAo|duia seysaul yusuleta ONULT tmory AQpaeo usstrq *Z-VII "Bly ysaw juawiaja [2x OZ (q) Ysa juaWa|a OZ X OZ (0) a ee ee so ; — T ] : i 4 I 2 rs) cC- 1 2 O=2n Qz=en {=!n =!n jUuawWaja saus09 4D juawiaja 49us09 Jo udlipuod AIDpuNog jo j1Dyaq UoljIpuod AuDpunog 40 j1DJaq -151- *saTtjord Ayyoojea aueldply sory Aqiaeo ueATIGg “*¢-y bly O'l- ysaw |2xOd ——— yseW O2x0Z -——- OOv= 24 Oe O'l- a ate | ysew {2x02 —— ysaw OZx OZ —-— OOl= 94 01 -152- ‘jusWala apou-% & Oy uoTWIpuos Arepunog rau103 snonutjuoostp jo uotzeurtxoiddy :moqf AqyiAeD uaATIGg *b-VWi B14 it -153- *ganoyuos AWTOWA0OA (p) pue {s1n0}U0d sins -saad (9) ‘souljuresrjs (q) ‘s10z20a AyToOTaA (ze) (YsOUT [ZX O2 ‘0 = 2u) MOTF AzTAEO U aaliq °*G-VIE ‘bid (P) (9) (9) (0) i a A a a ae a OO Ce eer 7 ot 64 Mm@ mma Lt MELLO COOL ITS NNN NE Lt lis fren SN NNN rit ts e7erer wT NSN NAAN rE seer, 2 22S NNNA\GV I phi fris. . OS NNNAN\\I Epi rss 22S pA NE [pi tsss----2 2 7p NNN [LN < — —— + = — at en ¥ -154~ *si1no0}u0d AW1O1I}10A (p) pue ‘sinojuOs 31ins -said (9) {‘sampuresrjs (q) ‘s10z3a0 AQTDOTaA (®) . ; s(ysaut [Z X OZ {001 = 2M) MOTT AZTAeD UAATIG own (P) (D) A A A ee ie 1 ot ft el Tn ae fe le a NO A a tL sss. \ RL AAESSEE PhAST? CIT SN NANAAYEE Reg; cS SANE Pyv ll fte say bayyee \ Sse rrr rT yyrye Noreen ee yy | YVAN aa Te == = _= eee ~155- *sanojuos AytoyW10Aa (p) pue ‘s1no0zUoD aiInB -sead (2) ‘sauy~ueerqys (q) ‘srozyoan AVDOTaA (B) (yseur [ZX OZ 00h = Pu) MOTs AQyAeO UaATIq Y-yII ‘bly (P) (9) (q) a {\¢ (0) or ggbrtrrrtttt ot A LA AH, \. A eg NN AAAs LAETOLI IIS SSAA i hE ool ITN NNNN hkl 7o sl TN NNNVA sereeeearas | aa ar | i RL ar A Av A | Rpt ae ar Ae ae \ WoT es i meme Ye TN ae mee ~ = -156- variation in velocity is employed for the 9-node element, and thus the amount of fluid "injected" into the cavity is different. Discrepancies noted between midplane velocity profiles (see [24]) may be attributed to the different boundary conditions, and are not indicative of the respective merits of the elements. The velocity field, streamlines, pressure contours, and velocity contours for Re = 0 employing the 20x 21 mesh are presented in Figure IIA-5. The corresponding plots for Re = 100 and Re = 400 using the same mesh are presented in Figure IIA-6 and Figure IIA-7, respectively. II-B. Entry Flow in a Channel In order to demonstrate the necessity of the “upwind" scheme, entry flow in a channel is studied. The problem statement and the cor- responding finite element mesh are shown in Figure IIB-1. We obtained the velocity field by employing the steady Navier-Stokes algorithm described above. Four different integration treatments of the advection terms are used and comparisons are presented in Figure IIB-2 for a Peclet number of 150. For details of the upwind scheme, consult Chapter 3 and Ref. [4]. As can be seen, the upwind schemes are superior to the Gauss schemes. For this high a Peclet number, there is little differ- ence between the optimal and fully upwind schemes. A much higher Peclet number (1.5x 107) was Studied. The optimal and full upwind schemes give identical results, whereas the Gauss schemes plot off scale (see Figure IIB-3). -157- *ysoul PUaWIsTa sqTUTT pueuotdiaosep walqoig :jeuueYys e ut MoT3 AtjUq *|-gT BI4 Yse-W (q) IX g 2 ¢ Z 0 I I I 40 3S¢° ey 40g’ JUBWI}LD]S LWA|GOsd (0) 9 > lx KAJQWWAS 4 g———— <X I=¢ | <—o-¢ . UID] q-°d G‘= e—____| Q=2n 04987 ig —____| | =!n OSl= 4 =8y *UlLa4 WOTJIIAPe ayy JO SyUSUZeaI} UOTeIZE4UT JUSsTaTIp Inojz 1ojz sj[Nsar jo uostieduroyg :[auueyo e ut mops Aatquy *2-gTT Bld ly A ly 9. GS 6 fF 2 Iai 9G » € 2 1 oO | Ly | | ho A PF SSNVD 2x2 Pe | ne SSNVS Ld! QGNIMdN T1NA Id | —— OQNIMdN IVWildO td! -158- *Laquimu jeToIeg YSsTY ye sauIayos putmdn [Nj pue yeumjdo 10j sz{nsey sjauueyo e ut MoTZ AtjUq ‘¢-gy B14 " 9s) << G2" ,Olx SI " 8d O'| -160- II-C. Couette Flow This problem actually is one-dimensional. The exact solution can be found in standard texts [25,26]. The problem description and mesh are shown in Figure IIC-1 and results are shown in Figure IIC-2. This is a simple problem in which the convec- tion term is identically zero. A boundary layer develops along the lower edge and diffuses upward, forming a steady, linear velocity pro- file as t increases. II-D. Dam-Reservoir Problem It is very important for engineers to know the hydrodynamic pressure along the inclined upstream face of a dam, due to earthquake motions. Explicit analytic formulas for calculating total horizontal, vertical, and normal loads are due to Chwang and Housner [27,28]. Approximate solutions and numerical solutions of the same problems have been given in [29,30]. The problem statement is shown in Figure IID-1. The dam is assumed rigid. The fluid is assumed to be incom~ pressible and inviscid. The initial conditions are quiescent and at t= 0* the dam face is set in constantly accelerating motion toward the reservoir. Data for the problems are given as follows: u = 90; o0=1; L= 2; H= 13 A= 10; y = 13 At = 0.025; and T = 0.1 (4 time steps). Meshes and results for the 60° dam and 90° dam are shown in Figure IID-2. Pressures are compared in Figure IID-3 with the exact, potential flow solution due to Chwang [27]. As can be seen, the re- sults are in good agreement. Also from Figure IID-2, the pressure contours are orthogonal to the streamline contours. It is expected from potential theory even though the full Navier-Stokes equation is *‘uolydizasap ulaTqoid pue ysaul JUsUaTa ayIUTY iMoTF ayenoy *}-OT 614 SJGON 11V lV O= én O< } | | O>4 6 }ein , -161- <a y= | yz. At = .0625 m_—— EXACT O FINITE ELEMENT FE Q 4 @ O 2 .29 l l l °% 2 4 6 8 Uy Fig. 1C-2. Couette flow: Comparison of finite element re- sults with exact solution. |.O ‘uoT}dt19saq ‘:watqoid r110oAr1asai-wed *)-am 614 ~163- (uoibas pazijas9sip) MOAIBSAY Q = 644, Quis 1) x= DI wog -164- *(---) sautjweaajs pue (—) sinojuo. ainssaid (fj) pue ‘s10zoaa AqpooTaA (3) ‘Ysa JUueteTA ayiUTA (P) *(,06 = 6) *(---) seuTjwest4s pue ¢ } sino}ues ainssaad (3) pue ‘s10jza0Aa Ayyaop2a (4) ‘Ysoaul yuatieyea ayury (2) (09 «x @) ‘warqoid q1oAtasai-weg *2-CTL 614 (4) (2) \ , 4 \ 4 \ ‘ \ \ SAN A A yen nwnwn SN mwes NNN NY -165- *sq[Nseat }OeKS pue yUeUtaTa aytuyy jo uosyieduiog :wayqoid 1toAresaz-uwieg "e-am BJ (Hod) /d =%9 o09=8@ OXI, — . {ua 06: juaWa|a ayUl4 oO O'| -166- employed here. II-E. Hamel Problem A problem statement is depicted in Figure IIE-1. The Hamel problem of convergent flow in a channel ("inflow problem") has been considered recently by several investigators (see [3,4,24,31]). The exact solution can be found in [25]. For the particular mesh employed, a radial velocity profile, in accord with the high Reynolds number approximation to the exact solution, is set at the outer radius (r= 4) having constant value of 1/4. The Reynolds number thus may be taken as Re = t/(6v). The circumferential velocity at r=4 is set to zero. The outflow boundary condition at r=1/4 is assumed traction free. At high Reynolds number this approximates the situation adequately. Results for Re = 500 and 5x 107 are presented in Figures II£-2 and ITE-3. As can be seen, the correlation with the exact solution ts very good. (Pressures in Figure IIE-2 and IIE-3 are reported at the element centers and are “unsmoothed.") II-F. Flow over a Step A problem statement is depicted in Figure IIF-1. We used this problem to demonstrate the necessity of "upwinding" when there 7s a "hard" upwind-facing boundary condition (see upwind section of Chapter 3) and too coarse a mesh is used. In Figure IIF-1 we present results of a calculation performed with 9-node Lagrange elements. A (product) Simpson's rule was used to construct the mass matrix, whereas Gauss- Legendre rules of order 3x3, 3x3, and 2x 2 were used on the convec- tion, u and A terms, respectively. Data employed were: p=1, p= 200; On DB ODO DO O -168- ? e e -—— Re = 500 l l ! O 2 4 6 1.0 68/11 (a) 8r ® —@ eo 6F ? © — Gre - ° x 9=4.3 oO. or l l © 2 3 , 4 (b) Fig.ILE-2. Hamel flow: Comparison of finite element (e) with exact (——) results at low Reynolds number. © (a) Velocity; and (b) pressure. -169- u, i ON BR OD} DW OT © e © o— Re = 5x10" a oS” r= 93 l | | i r’ O 2 4 6 8 1.0 60/1 (a) eo & ° © x au i i } J | 2 3 4 r Fig. IE-3. Hamel flow: Comparison of finite element (#) with exact ¢ ) results at high Reynolds number. ‘(a) Velocity; and (b) pressure. -170- “UOTJEULLOJUT YSou JUsUTaTS e}uly pue uoTydiaosap ulatqoig :deys e 1aa0 MoT “|-4qbI4 a ONIN habbo mbubhmbinbanbiuwhbb who wklnhuk hh LLL ye eh ww N aw Lk.d dé iil. | 4 juawayja ° a a (Uol}!pudd dv er oodk, 4 \\ 2pou -v <— }a|4NO ) y (lady oy osat = ~~ O=/4 A/Q=24 0 =2n Dal -In O<@neln — (UOHIPUCD ya}UL) | =n =| 7S 4 444 4L4Lbh 44h +hAhAAhbhAhShth+bhhbhAS,bthbhhbhhhSsbs+b+S;S;+EtEi*SsS+*s«‘AOis:}:ishiéised “ dp ee eee —~ _ — _ — — — _— a fe ea ee a ee —, —_ sO lO SO -— N N 7 ~ NON = ad ~ — — t \ we ~~ oe ees — — — — —_ “ ‘ . \ ~ —. — — — — — — — — ~e _ a - “ = ‘ ~ ~~ _— — — — — 7 = t=.7 (n=10) ee ee — - - _ — — — | — — —_— a ne ee — —_ —_ _— _—_ -— be Le es se a ed a ee ~ mee a ae ” |. 7 ra —_ ~ as = _— — _— ~— 7 ; / ~ \ N ~~ ees LL — — = — ” \ “ - 1 \ ~ — es — - —- = = - ee t= 1.4 (n=20) ee ee ee ee ~ - - - - - j- ee ee ee eed as — _ _ — j—- pee ge gn ee a —_ pe ee eR — ee -— _— — _ —_— 7 ' : - ~ ~ ~ “oat fea fe _ ’ . N ~ - ‘ \ ~N ~ — es — L- p— ee _ ~ - ~ ~ ~— —_ - . ~ — — —_ j— t=2.1 (n=30) Re a a see - - - Pe Le eee CCG _— —s eee eee ag — |-— be ee “ —_— TTT TTT re ng eg Te — i rr a oe o a Sa ee ee eel sees t— ™Y a N . ~ ~ — a a ee -_— — so — 7 , © ’ - ~ ~ ~ ~~ el — _ 7 . s ~ - . ~ ~ — — —~ + — os - - ~ ~~ — _ ~ . - —~ _— = Fig. 0F-2, t=2.8 (n=40) Flow over a step (Re = 200): Finite element results for 9-node elements with Gauss-Legendre integra- tion of convection terms. ee pee ee — —_— i }—. ee ee a OO a ee ee eee eee eC ee a 4 meee es ee es ee Foto t? “ODP TrTITIIItcte ee ' . ~ ~~ Oe oe -— _-_lcl _ _ ae ' . . ~ ~ ~ — — — — —_ }—- ee ~- . . ~ — ee t=.749 (n=10) — — — _ — i ee — — _ _ _— — _— _— = a ee esa Pe ee —=e ee ee as a re A a a ee —_— — ee 4 Se eee coed ee / - ~ ~~ ee ee ee eed a ae é/ “ ~ ~ ~ ~~ ~~ ~~ > ore a oa ee ed ¢ , - . ~ ~ ~ a es Oe t=1.50 (n=20) ee ee se se _ ~ - — L fe ee” _* "=e Fe ee es es — — — —_~ a eee ee ee — — rr a AO a ee +—~» ae 4 en a eed Phe Soe oe ee Re ae A ‘ . - ~ ~ — ~ sO COU ad _ va 1 . : . . ~ ~ ~ — — —_— t— t=2.25 (n= 30) a eee ~ ~- - - b a oe ee ad a a ee —~ | ~ ee oe ee Poe eee ee —_—— —_— ae / - = > — — ee ae A - - ~ ~ —~ — —~ eee tos _ — — - ca - ~ _~ _ _ — —~ a -~ - - - - - ~ = - bt t=3.00 (n=40) Re ee a —_ _ = - Le ee ee eee a pee a a a ee eee Oe ee ee ee ee en a ee ee 4 —_ ee ere ee tne ee iene ee ie ieee See Tl weed = _— —_— —* “— / - = =~ ~ — —_— Saenene Sates ee eee en Seed eee ee a ' . - - ~ ~ ~ ~ we we ee hb Pe ¢ ~ ~ ~ ~ _ = t=3.75 (n=50) : Fig. ILF-3. Flow over a step (Re = 200): Finite element results for 4enode elements with modified upwind treatment of convection terms. Lem a ee ee |. — — — a ee ess — — eet OO et ig eg a ee re a ae ee Mo at ee —~—s io rrr ee rr a ea 4 a ee aoe el / - ~ ~~ ~~ Ke Kerr rr OP eee “ / ~ ‘ ~~ ~ ~ ~ _ —— — fms eae , ‘ . ~ ~ meee eee ee eee - - ~ ~ —_ —- — ea t=.756 (n=10) ed —~_— — —_—_ andl ——- ee — — - —_ — cc — ~ —_— or eee” ease — — — — _ — eee ee pe ee ee ee es —_—-~-"— rr eee OO ee a ee ee ee ee _-—clc eo Ohm Lee —_ Wa —_ — — a a ee ee ~~ - - / a ae er oe a aes - ~ ~ ~ ~~ >= Or Tr fesse — —_— —_— ae 7 ~ ~ ~ — — — = a - - - - ~ = —_— ~— — La t=1.51 (n=20) a es ae ae ~ = i ee aad 1 a — —_— —_— ee i i in nn ey ey ee ey — rr at eee eee SS rr a a / — _— _— — a a ee fm — — —_ a / + - - ~ ~ ~~ — — — _— t— —_—- — “ / - - ~ ~ ~ ~ ~ — ~ — Peet ¢ . ~ ~ ~~ _ = — - - . -~ - - - - ~ _ = t= 2.26 (n=30) ee + ie i L- pean — ase ae eee ees OR i i i i i in ne a ee ee ee —_ te OO ee ee eee — ae / —_—_ — oe — _ OO —_ ae / - - - ad ~ _ ~ —~ oo Te en ae . - - ~ ~ ~ ~ ~~ KR — —_— _ _ - ¢ ’ . ~ _ b- ees ~ - - - - - - - - L od t= 3.02 (n=40) ee a a — i ~~ — —S ee a ae" Se” eee Li — — — — we i we ie ei a eR A Oe ON ee 7 — ee — — — — 7 — — ee ee a a en bd = _ _— _— “— / . - - - Sal ~ ~ ™ — ‘aed =e _— a / - - - ~ ~ ~ ~ a ee eee, : . “ L ee - - - ~ ~ ~ - - - - L t=3.77 (n=50) Fig. IIF-4, Flow over a step (Re = 107): Finite element results for 4-node elements with modified upwind treatment of convection terms, -174- A = 107; y = 1]; and At = 0.07. As is clearly visible, "wiggles" appear upstream of the step. Similar results are obtained for the 4-node elements employing Gauss-Legendre integration of the convective term. It is felt that this problem demonstrates the inappropriateness of Gauss- Legendre integration of the convected terms under the circumstances de- scribed. As can be seen from Figures IIF-3 and IIF-4, results for the 4-node elements which employ the modified upwind treatment of the con- vection term at Re = 200 and 10’, respectively, are smooth, and the upstream "wiggles" are removed in both cases. The data employed in these cases were: (Re = 200) p = 1; p = 200; A = 10°; y= 13 (Re= 10’) uw=13 p= 10!3, and y = 1. Time steps were selected adaptively. II-G. Axisymmetric Flow through a Sudden Enlargement The problem description and finite element mesh used are depicted in Figure IIG-?. The domain and mesh are split at Section iii for pic- torial purposes only. Macagno and Hung [32] have obtained both experi- mental and numerical results for this problem. A fixed time step of At=0.5 was employed, and y=0.75 was used for the transient algorithm. The dynamic viscosity, p, was set to 1 throughout. Thus Re=p and the penalty parameter was taken to be 10’. One point integration of the > term was employed. The calculations were performed in three se- quences and were compared to experimental and numerical results due to Macagno and Hung. In the first sequence, the initial conditions were quiescent and Re= 30. The sequence consisted of 60 time steps and a steady flow was achieved after approximately 30 steps. The flow was used as the initial condition for the second sequence in which Re= 60. This sequence consisted of 40 steps and a steady condition was attained *uoljdiz1osap ulsayqoig {jUsUlasie[ua ueppns & Ysno01Y3 MOTJ oTAJaUTUTASTXY *)- 91 6) 4 Ysa jUdW9]d ayiul4 (q) ; == EE PEE a Se a oe oe ee in wm TY JUdWA}DIS Wa|qOdd (D) 1! (+914NO) -2n-l “" O=ensl O=ensin | t i 2 4 (Asjausshs Jo sixo) Q=en=!y | at I? . 1H tt | ! O=en-=!n (O=¢n ‘5 1-| =In) t 4 a/\= ay ADO}aA yayul Id yo ae y_ cA (AsgamwAs jo sixD) O=e@nely | t | ~176- UW Uti, (a) ——— (b) +... _ >_> — THI, TUT ttttee TUTE thie Ltt Hittite. WHAM, WH ttttas. ttt. Wt. LUTE thre. TT ttttern. JUTE tunie. THU, TT tties, Lettie. THT ttttees. Mts, Hittites. Tuite. (tttiites. TH ttre. [ttteaas. Witte. Mitts. Itt. [Ittrs. [ftra,. Thine. Iie. ‘ Maa. Mi. Now] Wie, Thi, . Ihe. lhe. Meeees {b) vorticity con~ (d) streamlines; and (e) detail of streamlines in recirculation region. (e) Axisymmetric flow through a sudden enlargement (Re = 60): (d) (c) (a) Pressure contours; tours; (c) velocity vectors; Fig. IIG-2, — — = ‘UOISST UOTZE[NOAITIAT UL saUuTTUIeaI}s Jo [1eJep (a) pue ‘sauljuteerys (p) ‘sroyoadA Azio0Taa (9) fsino4 -uo0d AzioTIOA (q) ‘sinozuoo asinsseig (2) (007 = 24) yusulesie[ue Usppns e ysno1y} Mo]; olZeUIUIASTXY ‘*¢-OT BI 4 (0) ~178- after 20 steps. Results for this flow are presented in Figure IIG-2. With this flow as an initial condition, the final sequence in which Re = 200 was run for 90 steps. It took almost al] this time for the region just upstream of section iii to become steady. Results are pre- sented in Figure IIG-3, The overall flow patterns and length scales of the trapped annular eddies are in good agreement with [32]. To study flows at higher Re would require refinement and extension of the mesh downstream of section iii, as the recirculation regime tends to stretch out considerably with increasing Re. II-H. Viscous Flow about an Airfoil Viscous flow about an NACA0018 airfoil geometry is studied. The problem statement is shown in Figure IIH-1. The finite element mesh is shown in Figure IIH-2. We employed unit chord length (i.e., L) through- out. The dynamic viscosity, uu, was set to 1. The input inlet velocity is 1, therefore Re = o. Several runs were made with this mesh. The first, at low Reynolds number (i.e., 400) and quiescent initial conditions, shows well the diffusion of vorticity with the development of the boundary Tayer; see Figures IIH-3 and IIH-4. Accompanying pressure profiles are shown in Figure ITH-5. Steady, high Reynolds number results with quiescent initial conditions are presented in Figures IIH-6 and IIH-7. Obviously, from the velocity vector plots, this mesh is not adequate for such a high Reynolds number. The infinite-domain, potential flow, pressure profiles are presented in Fig. IIH-7 for comparison purposes. It is conjectured -179- *qUaUIaze}s Wal[qorg i]IOJATe ue JNOGe MOTF SNOOSTA “}-HT By t—_—_—_—162¢1| 1 a 10vE | ly <a O=en=\y e— O=en='4 A O=ene!n ¥ex (udIpuod (UOI}I DUDS }a]4NO) 4a|U!) 222" O=%4 a/=34 O=¢en o=4 l -In O=en=|] y -180- ‘YSoUL JUSUIaTS SPIUT | :[TOJATe ue ynoqe Mo]jJ snoostA “*g-yT'BI4 WVLAG + -->—~7 — — Pane = bn | ieee ee t= 3.0986 (n=280) Fig. IH-3. Viscous flow about an airfoil (Re = 400): Velocity vectors. ~182- f t=.2098 (n=20) t=.3299 (n=30) t=.4493 (n=40) t= 1.3293 (n=120) t= 3.0986 (n=280) Fig.ILH-4. Viscous flow about an airfoil (Re = 400): Vorticity contours. -183- ‘saainsseig :3(QQP = 2Y) [foyazTe ue ynoqe moj} snoost, ‘¢-HIL B14 -(O8Z=4) 60'S =} —— (OO2=¥) I2°2=4 ---- (OB =U) GQ =} verre OOv=9y ‘fuawala dasiul4 8'°O- vO- vO d/dz 8'0 c | -184- *sanojyuod ainssaid (2) pue ‘sautjuresays (q) ‘s10}00A AVPDOTIAA (®) *(901 = 94) [foyate ue ynoqe MoTZ snoost,A "g-HT by -185- 0.8 infinite-domain potential-flow solution —-— Re=10°, t= .5559 Finite element 0.6} an ee +=,5590 (7 (n-go) -198 ;- 2p/p Re=108, t = .5573 -0.8 Fig. IH-7, Viscous flow about an airfoil: Comparison of pres- sures at different Reynolds numbers. ~186- *sino}JUoO. 3ins “said (2) pue ‘sauqtyurearjzs (q) !s10jz00A AqTD0[9 A (8) 72°T = 338 (OL 93 OOP Worl OY Bursueyo Tazje paulezqo sz[nsay :[eosate "8-H B14 (9) -187- "901 °F 007 MOTT OY sutsueyo 1a4je peureyqo sainssaig :[Iojate ue ynoGe MOTT snoostA ° . 6-H 614 e 47'0 — =a TT thee ee, eit TS Te 42'0- ee. . bee Se eenientenl ah eed ae teeny wee et —_— ~ See tee Mech ee tee cae Sete ecm OY _ -” 2 = —— — —s ee —_— oe oe — —— ee (5901 9 OOD WO} 84 BuiBbubyod 43}4D) ( @oereavee juawea|a aflul4 2 O dz vO 9°0 -188- that enlarging the domain and refining the leading-edge region of the mesh would bring the results in closer agreement, although the pressure drop at the trailing edge is to be expected in a viscous computation. 6 calculation was made using the Re = 400 solution at An Re = 10 time 3.0986 as an initial condition. In this case the flow separated; sample results are shown in Figures IIH-8 and IIH-9. At later times (not shown) the separation point and recirculation regime moved down- stream. II-I. Axisymmetric Flow around a Sphere The problem description and mesh are shown in Figure ITI-1. Initial conditions were assumed quiescent and the "mean-incompressible" treatment of the A-term was employed. Runs at Re = 10 and 40 were made; sample results are shown in Figures I1 I-2and II I-3. After the flows became steady, comparisons were made with the infinite-domain, analytical results of Dennis and Walker [33]. As can be seen from Figure II I-3, vorticity (zc) and pressure are in good agreement. It is somewhat surprising that the results at Re = 40 are in such good agreement, as our finite-domain model seems hardly adequate for this high a Reynolds number. In Table II-I, pressure-drag (C,)s viscous-drag (C\). and total drag (Cp) coef- ficients are compared with the results of Dennis and Walker. {The coefficients are normalized as follows: drag/(mo(D/2)*.) The agreement is quite good overall, but better at Re = 10 than at Re = 40, as may be expected. ~189- ~™ N\ \ u)=l ug=O (inlet condition) \ \ \ \ ty 2up=O ty =Uo= 0 t=0 Ue20 Re=D/v 5.12D (outlet condition) aT ™ / plot region / H U)=U2=O0 —— 4.62D D “t | z i 4.62D (a) Problem statement ~ (b) Finite element mesh description and finite element mesh. Fig. 11-1, Axisymmetric flow around a sphere: Problem (c) Fig. 11-2. Axisymmetric flow around a sphere (Re = 10): (a) Velocity vectors; (b) pressure contours; and vorticity contours. ‘ ‘aanssaad (q) pue ‘AqypoW40,A (2) *sz[nsoir (——) Teo ATeue pue ( x ‘®) JuOWaTS o}TUTy jo uostzedwog :ersyds & punore mojyj olajauutAsTxy ‘“¢-] BIS (q) (D) (saaibap) g (saasbap) g O8!| OSI O¢ci O06 O9 OF OO O8| OSI O¢!l O06 O9 OF O | | 1 | | | | q | = + |- — qe ' — a — 1 -192- Table II-I Comparison of Drag Coefficients for Flow around a Sphere Re = 10 Re = 40 Present Dennis & Present Dennis & Coefficients Study Walker Study Walker Co 0.774 0.785 0.347 0.368 Cy 1.419 1.427 0.520 0.536 Cy 2.194 2.212 0.868 0.904 III. Remark on Slightly Compressible Formulation for Viscous Fluid Flow The slightly compressible formulation for viscous fluid flow (see Chapters 2 and 3) has. been programmed and studied. The dam-reservoir (IID), Hamel flow (IIE), and the free-surface wave generation flow prob- Tem (see Chapter 6) have been solved. The results obtained using this formulation are virtually identical with the penalty function formulation described above. So the sample examples will not be presented here. One of the main advantages of this formulation over the penalty function formulation is that when a static solution is reached, all accelerations and velocities will approach zero (due to the fact that it is a rate type constitutive equation, see Chapters 2 and 3). In the penalty formulation, a small nonzero velocity is produced to create a bulk viscous effect to support the loading (e.g., gravity). Hence, the slightly compressible approach is preferred and has been developed for the purpose of fluid-structure interaction problems which will be discussed in great detail in the next chapter. 10. -193- References T.J.R. Hughes, ‘Ae/Am 108 - Finite element method course notes," California Institute of Technology. T.J.R. Hughes, R. L. Taylor, and J. R. Levy, "A finite element method for incompressible viscous flows," in 2d Internat. Symp. on Finite Element Meth. in Flow Problems, S. Margherita Ligure , Italy, June 14-18, 1976. T.J.R. Hughes, R. L. Taylor, and J. F. Levy, “High Reynolds number, steady, incompressible flows by a finite element method," in Finite Elements in Fluids, vol. 3 (John Wiley, London, 1978). T.J.R. Hughes, W. K. Liu, and A. Brooks, "Review of finite element analysis of incompressible viscous flows by the penalty function formulation," J. Computational Phys. 30, 1-60 (1979). T.J.R. Hughes, “Implicit-explicit finite element techniques for sym- metric and nonsymmetric systems," Proc. Conf. on Nonlinear Problems in Mechanics, University of Swansea, Swansea, U.K., Sept. 2-5, 1980. M. A. Crisfield, "A faster modified Newton-Raphson iteration," Computer Meth. in Appl. Mech. and Eng. 20, 267-278 (1979). M. A. Crisfield, "Iterative solution procedures for linear and non- linear structural analysis," TRRL Laboratory Report 900, Dept. of the Environment, Dept. of Transport, 1979. K. C. Park, “An implicit, variable-step technique for fluid dynamics problems," in Preprints of 2d Internat. Symp. on Finite Element Methods in Flow Problems, S. Margherita Ligure, Italy, June 14-18, 1976, pp. 33-39. K. C. Park and P. G. Underwood, "A variable-step central difference method for structural dynamics analysis, Part I: Theoretical aspects," ASME Paper No. 79-PVP-120, presented at the Pressure Vessels and Piping Conference, San Francisco, CA, June 25-29, 1979. P, G. Underwood and K. C. Park, "A variable-step central difference method for structural dynamic analysis, Part II: Implementation and II. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. -194- performance evaluation," ASME Paper No. 79-PVP-121, presented at the Pressure Vessels and Piping Conference, San Francisco, CA, June 25-29, 1979. P. M. Gresho,, R. L. Lee, T. W. Stullich, and R. L. Sani, "Solution of the time-dependent Navier-Stokes equation via F.E.M.," 2nd Internat. Conf. on Finite Elements in Water Resources, London, England, July 10-14, 1978. S. L. Smith and C. A. Brebbia, Applied Math. Mod. 1, 234-266 (1977). T.d.R. Hughes, "Current trends in finite element research." Feature article to appear in Applied Mechanics Reviews. T.J.R. Hughes, "Recent developments in computer methods for struc- tural analysis," Nucl. Eng. and Design 57, No. 2, May 1980. T.0.R. Hughes and W. K. Liu, "Implicit-explicit finite elements in transient analysis: Stability theory," J. Appl. Mech. 45, 371-374 (1978). T.J.R. Hughes and W. K. Liu, "Implicit-explicit finite elements in transient analysis: Implementation and numerical examples," J. App]. Mech. 45, 375-378 (1978). T.J.R. Hughes, K. S: Pister, and R. L. Taylor, "Implicit-explicit finite elements in nonlinear transient analysis," Computer Meth. in Appl. Mech. and Eng. 17/18, 159-182 (1979). K. C. Park and C. A. Felippa, "Direct time integration methods in nonlinear structural dynamics," Computer Meth. in Appl. Mech. and Eng. 17/18, 277-313 (1979). E. Hinton,"Least squares analysis using finite elements," M.Sc. thesis, Civil Eng. Dept., University of Wales, Swansea, 1978. R. L. Lee, P. M. Gresho, and R. L. Sani, “Numerical smoothing tech- niques applied to some finite element solutions of the Navier- Stokes equation," 2d Internat. Conf. on Finite Elements in Water Resources, London, England, July 10-14, 1978. M. Bercovier and M. Engelman, J. Computational Phys., in press. 22. 23. 24, 25. 26. 27. 28. 29. 30. 3}. 32. 33. -195- 0. R. Bruggraf, J. Fluid Mech. 24, 113-151 (1966). D. K. Gartling, Finite element analysis of viscous incompressible flow, TICOM Report 74-8, University of Texas, Austin, December 1974. J. C. Heinrich, R. S. Marshall, and 0. C. Zienkiewicz, "Solution of Navier-Stokes equation by a penalty function finite element method," Report No. C/R/308-78, Dept. of Civil Eng., University of Wales, Swansea, 1978. G. K. Batchelor, An Introduction to Fluid Mechanics (Cambridge University Press, Cambridge, England, 1970). H. Schlichting, Boundary-Layer Theory (McGraw-Hill, 1968), 6th edition. A. T. Chwang,"Hydrodynamic pressures on sloping dams during earth- quakes. Part 2. Exact theory." J. Fluid Mech. 87, 343-348 (1978). A. T. Chwang and G. W. Housner, "Hydrodynamic pressures on sloping dams during earthquakes. Part 1. Momentum method." J. Fluid Mech. 87, 335-342 (1978). T. Von Karman, Trans. Amer. Soc. Civil Eng. 98, 418-433 (1933). S. K. Sharan and G.M.L. Gladwell, "A method of analyzing dam- reservoir-foundation interaction problems," Proc. Symp. on Appli- cations of Computer Methods in Engineering, Vol. II, pp. 1417- 1487, University of Southern California, August 23-26, 1977. D. K. Gartling, R. E. Nickell, and R. Z. Tanner, Internat. J. Numerical Meth. in Eng. 11, 1155-1174 (1977). E. 0. Macagno and T. K. Hung, J. Fluid Mech. 28, 43-67 (1967). S.C.R. Dennis and J.D.A. Walker, J. Fluid Mech. 48, 771-789 (1971). -196- Chapter 6 APPLICATION TO FLUID-STRUCTURE INTERACTION I. Introduction In Chapter 4 a general three-dimensional finite element formula- tion was described for quasistatic and dynamic nonlinear shell analysis. This shell theory can easily be degenerated to an analogous two- dimensional formulation which includes several special cases of prac- tical interest. Numerical examples indicate the good behavior of the elements studied. In Chapter 3, the penalty function/finite element formulations of the Navier-Stokes equations was described. Numerical examples indicate the effectiveness of this scheme. The slightly compressible formulation for viscous fluid flow has been programmed and studied. As has been mentioned, one of the main advantages of this formulation over the penaity function formulation is that when a static solution is reached, all accelerations and velocities will approach zero. Hence, this slightly compressible approach will be used for the purpose of fluid- structure interactions. As a practical application, the proposed finite element proce- dures are employed for the dynamic, three-dimensional, nonlinear, inelastic response of ground-supported, cylindrical liquid storage tanks which accounts for fluid-structure interaction and free-surface sloshing. Early studies relevant to this topic were performed by Westergaard [1]. Hoskins and Jacobsen [2], Jacobsen [3], Werner and -197- Sundquist [4], and Jacobsen and Ayre [5]. They assumed the tank to be rigid and the hydrodynamic wall pressures are then determined. Housner [6] developed an analytical method for determining hydrodynamic wall pressure under the same assumption that the tank is rigid. A mechanical equivalent model is developed and is widely used in current earthquake design practice [7]. Within recent years, studies performed by Veletsos [8], Yang [9], and Veletsos and Yang [10] have shown, however, that the rigid- tank assumption may lead to a significant under-estimation of the magnitude of the resulting forces. The first use of a digital computer in analyzing this problem with the aid of finite element methods was completed in 1969 by Edwards [11]. This investigation treated the coupled interaction between the elastic wall of the tank and the con- tained liquid. Recent developments along the same line were performed by Wu et al. [12] and Shaaban and Nash [13]. Haroun [14] recently conducted a theoretical and experimental] investigation of the dynamic behavior of cylindrical liquid storage tanks to seek possible improvements in the design of such tanks to resist earthquakes. In this study, natural frequencies of vibration of the associated mode shapes are found through the use of a discreti- zation scheme in which the elastic shell is modeled by finite elements and the fluid region is treated as a continuum by boundary solution techniques. However, the methods cited are based on smal? motion linear response, and do not take satisfactory account of nonlinear effects observed in the field under the action of strong earthquakes, -198- such as lift-off of tanks from foundations, finite-amplitude sloshing, and nonlinear and inelastic tank response. Obviously a more advanced approach is called for. The whole subject is indeed very complicated and it is not the purpose of this thesis to perform a nonlinear tank study, yet the author hopes to have developed a tool which can be used eventually for that purpose and others. Before we can perform a fluid-structure interaction calculation, a contact/sliding element is required to handie the siiding assumption of the fluid on the tank wall. We begin with the discussion of such a sliding element. A free-surface wave-propagation problem is then used to determine the effectiveness of the proposed mixed Lagrangian-Eulerian scheme. The classical static axisymmetrical buckling of elastic cylind- rical shells subjected to three types of boundary conditions is then studied, followed by a buckling analysis of a tilted cylindrical liquid storage tank subjected to internal hydrostatic loading. Finally, a dyn- amic analysis of a liquid-filled cylindrical tank is performed to demon- strate the finite element fluid-structure interaction procedures. II. Contact/Stiding Element Contact elements may be used to impose inequality constraints between nodes. Either perfect friction (i.e., “stick") or frictionless (i.e., "“slide") conditions may be achieved. A contact element is defined by three nodes and a spring constant, or "penalty parameter," k, for our purposes. The connection from node A to node B defines the "slide line" direction; whereas, node C is the contact node (see Fig. II-1). The normalized projected distance of node C to node A is denoted by a € [0,1] and is given by: ~199- are a\\ Fig O-1I A contact element is defined by three nodes -200- > a = AC AB/ [ABI (11.1) > Here XY denotes a vector connecting point X to point Y and "+*"is the dot product. The normal vector 2 from node C to AB is defined by: rh > g =" (AB x AC) x AB (11.2) ~ a+ 2/ || (11.3) where "x" is the cross product and "+" means "is replaced by." The local contact stiffness matrix Ke is: (l-a)* a(l-a) -(1-a) g 2 k sym of -O (11.4) aw i ] The rows and columns of the contact stiffness matrix are arranged in such a way that the first row or column corresponds to node A, while the 2nd row (or column) and 3rd row (or column) correspond to nodes B and C, respectively. Before the global assembly procedure, the local contact stiffness matrix is rotated to the global form Ke (for the sliding case) via: Ke TK? tT (11.5) TT | £ 013 943 T=f0 yg oO (11.6) T=/0)3 & 03 3 0., 0 Q p13 S13 -201- 0 [fo oO 0] (11.7) 213 - For the sticking case KS is F(1-a)* 0 0 a(l-a) 0 0 -(I-c) oO 0 (1-a)* 0 0 a(l-a) 0 0 -(I-a) 0 (I-a)* 0 0 a(l-a) 0 0 (I-a) of =O 0 +a 0 0 Ke = a ~=s0 0 -« O . SYMMETRIC oe =O 0 -o 1 0 0 LL 1 J (11.8) III. Numerical Examples III-A. Free-Surface Wave-Propagation Problem The problem statement and finite element meshes are shown in Figure IIA-1. We generate a wave numerically by prescribing a displace- ment-time history at the left-hand boundary of the domain. This pre- scribed function is designed in such a way that it is consistent with the experimental studies performed at Keck Laboratory of the California Institute of Technology. Its function is to produce a solitary wave-- one in which the nonlinear and dispersive effects are balanced so that the wave propagates without distortion. Data used in the analysis are given as follows: g = 13 L = 949.095; D = 10; H = 0.86, h = 5.895; and At = h/c = 1.7888. Meshes 1 and 2 consist of 160 and 320, 4-node elements, respectively. The flow was asSumed inviscid, p = 1 and A = 2.6089x 107, Lt G00, Kost)? HE] + tanh (SET -4 4)]; c=./gD(ItH/D); cof oe | | | —— REFERENCE DOMAIN 2, ! ——— DEFORMED DOMAIN 2, A“ je ———» Po N “x L (a) Problem statement and schematic of domains rr L toe toe 7 = : (b) Finite element mesh | (c) Finite element mesh 2 m | i m Ke L- Lagrangian degree-of-freedom E- Eulerian degree-of- freedom Fig. If A-| Problem statement and finite-element meshes for free-surface flow -203- The finite element calculations employing a one-dimensional depth averaged theory [15] are compared with calculations for mesh 1, one element through the depth. The results are shown in Figure IIIA-2. The figure shows the results agree remarkably well, considering the two approaches are quite different. Notice that the rate of growth of the trough is a little greater with the present scheme than it is for the one-dimensional theory. Even though there is good agreement between the two theories, experimental results are not consistent with the previously obtained numerical results. The differences are: (1) the experimental waveform propagates somewhat faster and exhibits less amplitude decay; and (2) dispersion emanating from the back of the wave is not seen experimentally. The reason for the discrepancy between the numerical and experi- mental results is the variation of the velocity through the depth is not linear. Both the one-dimensional depth average theory and one element through the depth, two-dimensional theory invoke this assump- tion. As can be seen in Figure IIIA-3, with two elements through the depth, the trough has been reduced and the relative wave height is essentially constant with propagation. The difference in having two instead of one element in depth is that the distribution of velocity is no longer constrained to be linear with depth. It is also noted experimentally that "noise" superposed upon a solitary wave tends to propagate at a slower velocity than the main pulse. A numerical simulation of this phenomenon is shown in Figure IITA-~4. The noise is generated by setting a no-slip boundary condition typical Us (x,,D,nAt)/D —--- one-dimensional, — depth- averaged 7 finite element code two-dimensional Lagrangian- Eulerian formulation | (mesh |-one element_ through depth) L ! ! 40 50 60 100 x,/D ‘Fig. IT A-2 Comparison of one-dimensional depth-averaged theory and two-dimensional Lagrangian - Eulerian formulation _ f.l =i" a So; £ a -— < = _ a a ~ =|0h—~ pnamey 4 = 4 i 4 - — a _ > 7 <J _ > | = ay x a Ss _ cof a woneneoe mesh | (one element -— through depth) 7 mesh 2 (two elements 7 through depth) 4 I l a J uo 60 80 100 x,/D Fig. II A-3 Comparison of results for meshes | and 2 -206- oIb Older sy 4 a — OF amet A 4 Ss a <] a c | 4 fa} _ et Z = J <s - wnneoa-- mesh | (one element _ through depth) 4 mesh 2 (two elements 7 through depth) 4 ———— mesh 2 with superposed- “noise a n=20 _ t l | \ Ld 0 20 40 60 80 100 x,/D Fig. IZ A-4 Comparison of results with, and without, superposed ay . u“ noise , -207- along the left-hand edge of the domain, in contrast to previous calcula- tions in which a slip condition was employed. Note how the solitary wave emerges, unaffected by the initial perturbations. III-B. Static Axisymmetric Buckling of Elastic Cylindrical Shells Three types of boundary conditions were employed in this study: a) bottom fixed--top free; b) bottom fixed--top fixed with respect to horizontal displacement and rotation; and c) bottom simply supported--top fixed with respect to horizontal dis- placement. Load versus displacement data and deformed profiles are shown in Figure IIIB-1 up to the buckling loads at which time the calculations were ter- minated. Comparison is made with analytical results based upon eigen- value analysis (see e.g., [16]). I1I-C. Cylindrical Shell Containing a Fluid An analysis was performed of a tilted cylindrical shell sub- jected to internal hydrostatic loading. The problem statement is shown in Figure IIIC-1, along with the stress results prior to any significant nonlinear phenomena. Data employed in this analysis are given as fol- lows: E = 7.35x10° 1b/ine (Young's modulus) v= 0.3 (Poisson's ratio) D = 16.0 in L = 20.0 in R= 4,0 in h = 0.01 in a = 30.0° -208- “S|J8YS [DI4PUAD Diysoja jo BHuijyonqg dWspewwAsixo d4oOIS ‘am bi4 1s) S| Vv O°Ob S| O'lz 4) ’ . G*0d $ . . 0'0Z a , , ony | Pp ome) ge O'S] 9 Gd 2 omel| q O'S | O'S D (| <Ol) pool |! youd (G1 .01) Poo] alyoud (Ulz_O]) }UsWAadD]dsip |DdI,4aA Id gl Sl ral 6 9 ¢ ie) T T T T T T fe) SLUSWS|9 SPOU-~F7 Ajuamy B nian UdlyNnjos 70! . pabsaauod isp ® 4 IE | 402 =>. 4 Ky *{DOISSDIO y —- a LOL ec’ =a zui/qi ,OIxX¢e =4 | A ul So°O =Y 4dOv Je Eg wOl = ee ge ee ee ulg=4 ] ! j | (Qi,O1) Poo] 1040} -209- *“awiBes sDaul] suf Ul S4|NSA1 SSBatS “Bulpoo} 91;DySO4pAY |DUJa}U! Of pajoalqns jjays JOospuljAS paxil “b- OMT GI4 ‘ lo \-O1 x (NI) ®X JONWISIO WIx a= iy &x SA SSBULS TVIXY S 9" rad 80 v0 (exe) fe) ‘o 1 T T T or b ro) x g oF de om w wo a -4 a n m Lim vB on IN wm oO. NvaN —— ° 3 ZISLNO +++ 3 a aaisNi —*—*- of NYSA wn to X AGISLNO ee Ho x pan 3disNi ~*—e— ° © st S) ¢ He aes xX SA S$S3YHLS dOOH 9" ra 80 v0 oo? % AV s 00 |-O1X (NI) &X JONVLSIO WIXV t 3 2 SYNLVTONSWON ONY SNOILIGNOD SNIGVO) Kot BS x > ' T N 10 1-01 X (NI) ®X 39NVLSIC TIX 4 ot 2 8'0 v'0 “gol od i, a ley —Te T vr * *T w o> | = . in 2 x \ fo} x f 7 = Ho wi a= yy &X SA INIWOW TWIXV o usd O=¢ jo &y -210- y = 0.0361 1b/in® (density of water) The effective specific gravity, Vorp> Was defined to be Yore * Y9- The value of g used in obtaining the stresses shown in Figure IIIC-1 was 1.8. Comparison with a membrane theory solution [17] confirmed the cor- rectness of the numerical results outside of the boundary layer at the fixed end. Buckling calculations were performed by incrementing g. Results are shown in FiguresIIIC-2 and IIIC-3. As may be seen, the initial bifurcation occurred in attempting to increment g beyond 1.8. This corresponds to a maximum axial membrane stress of 547.53 psi at x = L, o = 7, which is approximately half the classical value. Although this result was not anticipated, it bears a great deal of similarity to results obtained by Argyris and Dunne [18] in their penetrating study of a compressed cylindrical panel. They too found an initial bifurcation of approximately half classical. When they pursued the analysis into the post-buckled regime, a drop in load of only 3% was noted, at which time the load was able to be increased again to a value in the vicinity of classical. At this point a second, and more significant bifurcation ensued, and this was followed by several more bifurcations. An impor- tant point in considering the results of Argyris and Dunne was that a preliminary eigenvalue analysis, employing the same mesh used in the finite deformation case, yielded a buckling load approximating the classical value on the high side. This lends credence to the analysis, but a full understanding of the computed initial bifurcation in this case, and the present case,does not yet seem to exist. The results shown were obtained using $1 elements. Similar results were also obtained -211- FH yall } 1 WAL / 7; 7 [ aaa [ ; | ma Noy Tr he ies bh bh = = ~~ A | he i a = ix ‘a ~ 1 rm a | hy ~ ht = = a Na Hy = = wtp IM At tp} seme —_ ie i os — th — ots t-—t g= 0.00 fe) gzl. rf { an 44 L\ \ fe rH VA \ SS i wa fi +- / Ly | i ee \ ine (| MY line R 4 rt Ni ht | i hina — ssSS= = === g=1,.80 g=1.82 g=1.90 (START FROM g=1.80) | (START FROM g=1.80) Fig. IIC-2, Deformed configurations for the tilted ‘cylindrical shell. -212- g= 0.60 —— g2!.90 (START FROM g=1.80} SSS 00 g=0 g g 7182 (START FROM q=1.80) Fig. IIC-3. Detail of lower regions for the tilted cylindrical shell. -213- with heterosis elements. III-D. Dynamic Analysis of a Liguid Filled Cylindrical Tank A dynamic analysis was performed of a tiquid filled cylindrical tank subjected to a periodic saw-tooth ground acceleration excitation. The problem statement is shown in Figure IIID-1, along with displace- ment, velocity, and acceleration time histories at three different positions (node A, B, and C). The tank is subjected to the mentioned excitation in the x-direction only. Therefore all the bottom nodes of the tank are subjected to the same excitation with the other two direc- tions (y and z) fixed. A typical input ground excitation (acceleration, velocity, and displacement) is given for node A. Due to symmetry, only half the tank is modeled by 36 S1 shell elements, 80 uniform-reduced in- tegration 8-node fluid elements, and 30 sliding elements. The height of the tank is 864 ins (i.e., distance between A and C). The mean dia- meter of the tank is 577 in. (i.e., distance between B and C). The tank thickness is 1 inch. It is filled with water up to 720 in. The material properties of the tank and water are: Tank (steel): E = 3.0x10" psi v= 0.3 0 = 7.324x 107" Ib sec“/in’ oy = 3.0x 10" psi (yield stress) B = 0.0 (plastic modulus) Water: B = 6.991x10!* psi (bulk modulus) n= 0.0 (dynamic viscosity) -214- ee ee {Mt 01} i ¥ ONY 9 JOON N3SM138 LNIW30V SIO JAILY TI Co a ee oa ay) 2 300N 30 LNGINSOUTdSIC eee 8 = or NL 7 _ < i , eee ~ 7 ‘le : ce a nm (oma ony 3 KON 40 ALIDOTIA (pers on 3 3008 40 NOLLWER TOW a ee Te eo (4 on ¥ ONY @ 300N N33K138 “INSH2DV 12510 SAIL 38 Td cat eee ORT Raa wae g bi bs Hy L oe Cla cl a ed (uy) 8 300N 30 LNB TSI smerny. a RAM AY LT FINE ELEMENT MESH ¥ 300N 4 ALISOIIA Fig. IIID-1. ~215- Fig. I D-2. -216- maaan it OES ll Fig. ID-3. -217- 9.349 x 107" 1b sec/in® 2 aD u 3.864 x 10° in/sec (gravity) oO ili The spring constant, k, of the contact elements is taken to be 4.44 x i9!6 psi. The Newmark parameters are: 8=0.5625,y=1.0, andAt=0.02 sec. A linear implicit-nonlinear explicit operator splitting algorithm with one iteraction is employed. A small deformation analysis option is used. The maximum computed stress occurring in the tank is less than oy even though the peak ground excitation is 1.5 g. The relative free surface motion as well as the relative tank displacement is within the small motion region. Therefore, no further analysis of this problem was performed. The computed acceleration, velocity, and displacement time histories of nodes B and C as well as the relative displacement between nodes B and A and between nodes C and A are Shown in Figure IIID-1. It is expected that the fundamental mode is dominant for this kind of excitation; also nodes B and C are expected to respond in "phase." This phenomenon can be seen fron the absolute displacement- time histories of nodes B and C in the figure. The deformed shape of the tank and the free surface motion (with the input translational motion subtracted) are shown in Figure IIID-2 at a time interval of 0.2 seconds. The magnification factor of these plots is 7.81. Frequency responses higher than the fundamental slosh- ing frequency were observed. The sloshing period is approximately 6 seconds. Since the free surface motion is much larger than the tank deformation, plots of the tank response (with the input translational motion subtracted) are shown in Figure IIID-3, magnified by a factor of 76.31. 10. 11. 12. -218- References H. M. Westergaard, "Water pressure on dams during earthquakes,” Trans. ASCE 98 (1933). L. M. Hoskins and L. S$. Jacobsen, "Water pressure in a tank caused by a simulated earthquake," Bull. SSA 24 (1934). L. S. Jacobsen, “Impulsive hydrodynamics of fluid inside a cylind- rical tank and of a fluid surrounding a cylindrical pier," Bull. SSA 39 (1949). P. W. Werner and K. J. Sundquist, "On hydrodynamic earthquake ef- fects," Trans. AM Geophys. Union 30 (1949). L. S. Jacobsen and R. S. Ayre, "Hydrodynamic experiments with rigid cylindrical tank subjected to transient motions," Bull. SSA 41 (1951). G. W. Housner, "Dynamic pressures on accelerated fluid containers," Bull. SSA 47 (1957). G. W. Housner, "The dynamic behavior of water tanks,” Bull. SSA 53, (1963). A. S. Veletsos, "Seismic effects in flexible liquid storage tank," Proc. 5th World Conf. on Earthquake Eng. 1, 1974. J. Y. Yang, "Dynamic behavior of fluid-tank systems,” Ph.D. thesis, Rice University, Houston, Texas, 1976. A. S. Veletsos and J. Y. Yang, “Dynamics of fixed-base liquid- storage tanks," Proc. U.S.-Japan Seminar on Earthquake Eng. Research with Emphasis on Lifeline Systems, Tokyo, Nov. 1976, pp. 317-341. N. Edwards, "A procedure for the dynamic analysis of thin walled cylindrical liquid storage tanks," Ph.D, thesis, U. of Michigan, Ann Arbor, Mich., 1969. C. I. Wu, T. Mouzakis, W. A. Nash, and J. M. Colonell, “Natural frequencies of cylindrical liquid storage containers,” Dept. of Civil Eng., University of Massachusetts, June 1975. 13. 14. 15. 16. 7. 18. ~219- S. H. Shaaban and W. A. Nash, "Response of an empty cylindrical ground supported liquid storage tank to base excitation," Dept. of Civil Eng., University of Massachusetts, August 1975. M. A. Haroun, "Dynamic analysis of liquid storage tank," Ph.D thesis, California Institute of Technology, Dec. 1979. . D. G.-Goring, "Tsunamis--The propagation of long waves onto a shelf," Ph.D. thesis, Report No. KH-R-38, California Institute of Technology, Nov. 1978. S. P. Timoshenko and J. M. Gere, Theory of Elastic Stability (McGraw-Hill Book Co., New York, 1961). C. Shih and C. D. Babcock, Jr., Private communication. J. H. Argyris and P. C. Dunne, "Post-buckling, finite element analysis of circular cylinders under end load," Report 224, Institut fiir Statik und Dynamik der Luft und Raumfahrtkonstruk- tionen, University of Stuttgart, Germany, 1977. -220- Chapter 7 SUMMARY AND SUGGESTIONS FOR FURTHER DEVELOPMENT In this thesis, finite element procedures for fluid-structure interaction problems have been developed. Different aspects of non- linear methodologies have been studied. A working finite element method computer code [15], which will be discussed in a separate report, has been developed which may be applied to the dynamic, three-dimensional, nonlinear, inelastic analysis of ground-supported, cylindrical liquid storage tanks subjected to strong shaking, and various other fluid- structure phenomena. This is a significant step toward the development of a clean general purpose, modularized, fluid-structure interaction, finite element computer program that engineers will find helpful. Future fluid-structure analysis developments will need to include improved transient algorithms, better fluid-structure interfac- ing techniques, and contact-impact techniques. Automatic time stepping strategies based upon accuracy consid~ erations [1,2], subcycling techniques [3-5], and effective iterative equation solvers [6-9] will be generalized and eventually included in this fluid-structure computer program. The presently developed, slightly compressible fluid formulation is a step towards the inclusion of a compressible formulation, although further studies are clearly required. The mixed Lagrangian-Eulerian method seems to treat fluid- structure interaction problems with large motions of the structure naturally and economically. However, for large-deformation bubble ~221- dynamics calculations [10] that occur, for example, in light water reactor systems, “continuous rezoning" of meshes, together with mixed Lagrangian-Eulerian methods are required. As rezoning is expensive and difficult to do adaptively, other techniques are called for. Hirt and Nichols [11] proposed a volume-of-fluid finite difference method (VOF) in which the boundary of the free surface is not specified directly, but instead is defined by the fraction of fluid in each cell. The fraction of fluid is treated as a cell variable and governed by a transport equation. This technique can permit the solution of very complicated free-surface problems (e.g., bubble dynamics, etc.) , It is not yet clear how to develop VOF finite element methods due to the more general topology of finite element meshes compared with finite differ- ences. A combination of VOF and mixed Lagrangian-Eulerian methods may be an effective compromise. The contact/sliding element proposed here for fluid-structure interaction problems cannot be used directly for tanks lifting off from their foundations, and the resulting large-amplitude sloshing causing roof damage. This contact element, together with the contact-impact technique developed by Hughes et al. [12,13] and Hallquist [14] can serve the purpose and would be a worthwhile generalization. Finally, a word of caution is in order regarding the ability of the present capabilities to solve general problems of fluid and struc- tural response. Although a degree of confidence has been established in certain realms of phenomena, it should be emphasized that the present developments are incapable of solving many problems of -222- physical interest and application of these techniques in such circum- stances is speculative at best. Areas which are particularly difficult and presently beyond the developed capabilities are (among others) high Reynolds' number flows and complex shell buckling. Much numerical research still needs to be done in these areas. 10. -223- References K. C. Park and P. G. Underwood, "A variable-step central difference method for structural dynamic analysis, Part I: Theoretical aspects," ASME Paper No. 79-PVP-120, presented at the Pressure Vessels and Piping Conference, San Francisco, CA, June 25-29, 1979. P..G. Underwood and K. C. Park, "A variable-step central difference method for structural dynamics analysis, Part II: Implementation and performance evaluation," ASME Paper No. 79-PVP-121, presented at the Pressure Vessels and Piping Conference, San Francisco, CA, June 25-29, 1979. T. Belytschko, H. J. Yen, and R. Mullen, “Mixed methods for time integration," Computer Methods in Appl. Mech. and Eng. 17/18, 259- 275 (1979). J. P. Wright, "Mixed time integration schemes," Computers and Structures 10, 235-238 (1979). J. Isenberg, D. K. Vaughan, I. Sandler, "Nonlinear soil-structure interaction," EPRI NP-945, Project 810-2, Final Report, Dec. 1978. M. A. Crisfield, "A faster modified Newton-Raphson iteration," Computer Meth. in Appl. Mech. and Eng. 20, 267-278 (1979). M. A. Crisfield, “Iterative solution procedure for linear and non- linear structural analysis," TRRL Laboratory Report 900, Dept. of Transportation, 1979. T.J.R. Hughes and J. M. Winget, "Symmetric and nonsymmetric equa- tion solvers," to appear. T.J.R. Hughes, "Implicit-explicit finite element techniques for symmetric and nonsymmetric systems," Proc. Conf. on Nonlinear Problems in Mechanics, University of Swansea, Swansea, U.K., Sept. 2-5, 1980. J. Donea, P. Fasoli-Stella, S. Giuliani, J. P. Halleux, and A. V. Jones, "An arbitrary Lagrangian-Eulerian finite element procedure for transient dynamic fluid-structure interaction problems,” paper No. B1I/3, Trans. 5th Internat. Conf. on Structural Mechanics i. 12. 13. 14, 15. “2? 4- in Reactor Technology, Berlin, Germany, Aug. 13-17, 1979. C. W. Hirt and B. D. Nichols, "Eulerian method with free boundaries, in Preprints of the Ist Internat. Seminar on Fluid Structure Interactions in LWR Systems, International Congress Center, Berlin, (West) Germany, Aug. 20-21, 1979 (T. Belytschko, ed.), pp. 29-42. T.J.R. Hughes, R. L. Taylor, J. L. Sackman, A. Currier, W. Kanoknukulchai. "A finite element method for a class of contact- impact problems," Computer Meth. in Appl. Mech. and Eng. 8, 249- 276 (1976). T.d.R. Hughes, W. Kanoknukulchai and R. L. Taylor, "A finite ele- ment method for large displacement contact and impact problems," in Formulation and Computational Algorithms in Finite Element Analysis (M.1.T. Press, Cambridge, 1977). J. 0. Hallquist, "A numerical treatment of sliding interfaces and impact," Computational Techniques for Interface Problems, AMD 30, ASME, 1978. W. K. Liu, “FLUSTR, A fluid-structure interaction finite element analysis computer program," User's manual, in preparation.