PROGRAMMING THE FINITE ELEMENT METHOD 2ND EDITIONPROGRAMMING THE FINITE ELEMENT METHOD SECOND EDITION LM. Smith and D.V. Griffiths University of Manchester, UK JOHN WILEY & SONS Chichester . New York . Brisbane . Toronto . SingaporeCopyright © 1982, 1988 by John Wiley & Sons Ltd. All rights reserved [No part ofthis book may be reproduced by any means, for transmitted, oF translated into a machine language without the written permission of the publisher. Library of Congress Cataloging-in-Publication Data: miming the finite element method, Includes bibliographies and indexes. L. Finite element method Computer programs, 2 Engineering - Data processing 3. Soil mechanics—Data processing, 1. Griffiths, D.V. Ml Title TA347.FSS64° 1988 620.001'S15353 87-8123 ISBN 0 471 915352 1 (cloth) . ISBN 0 471 91353 X (paper) British Library Cataloguing in Publication Data: Smith, 1. M. (Ian Moffat) Programming the finite element method. —2nd ed 1. Finite element method—Data processing [Title I. Grilfths, D. V. 515353 TAMI.PS ISBN 0 471 91552 1 (cloth) ISBN 0 471 91353 X (paper), ‘Typeset at Thomson Press (India) Limited, New Delhi Printed and bound in Great Britain by Anchor Brendon Ltd, Tiptree, EssexContents Preface to Second Edition. Chapter 1 Preliminaries: Computer Strategies 10 Introduction : oo 1.1 Computer Hardware... 12. Store Management . rn 13. Vector and Array Processors 5... 14 Software . : 15 Portability of Subroutine’ or ‘Procedures’. 16 Libraries. bp oc omeos 17 Structured Programming | |). ) 1. 18 Other Languages. 2 19 Conclusions ee 1.10 References Chapter 2 Spatial Discretisation by Finite Elements 20 Introduction 21 Rod Element 22. Rod Inertia Matrix ee 23 The Eigenvalue Equation. 2). 24 Beam Element oo 25. Beam Inertia Matrix. ee 26 Beam with an Axial Force... . . 27. Beam on Blastic Foundation 28 General Remarks on the Discretisation Process. 2.9 Alternative Derivation of Element Stiffness . 2.10 Two-dimensional Elements: Plane Stress and Strain’ 211 Energy Approach : 2.12 Plane Element Inertia Matrix 2.13 Axisymmetric Stress and Strain 2.14 Three-dimensional Stress and Strain . 2.15 Plate Bending Element 2.16 Summary of Element Equations for Solids . 2.17 Flow of Fluids: Navier-Stokes Equations. 2.18 Simplified Flow Equations 2.19 Simplified Fluid Flow: Steady State2.20 Simplified Fluid Flow: Transient State . 2.21 Simplified Fluid Flow with Advection 222 Further Coupled Equations: Biot Consolidation . 2.23 Conclusions. 2.24 References Chapter 3 Programming Finite Element Computations 3.0 Introduction... 3.1 Local Coordinates for Quadrilateral Elements. 3.2. Numerical Integration for Quadrilaterals 33. Local Coordinates for Triangular Elements 3.44 Numerical Integration for Triangles 3.5 Element Assembly 3.6 Incorporation of Boundary Conditions . 3.7. Programming Using Building Blocks 38 Black Box Routines 3.9. Special Purpose Routines 3.10 Solution of Equilibrium Equations. 3.11 Evaluation of Eigenvalues and Figenvectors 3.12 Solution of First Order Time-dependent Problems . 313 Solution of Coupled Navier-Stokes Problems: 3.14 Solution of Coupled Transient Problems 3115 Solution of Second Order Time-dependent Problems 3.16 Conclusions 3.17 References Chapter 4. Static Equi 40. Introduction ibrium of Structures + Equilibrium of Uniform Beams Program 40 : Program 4.1 : Equilibrium of Stepped Beams Incorporating Prescribed Displacement Boundary Conditions Program 42 : Beam on Elastic Foundation with Numerically Integrated Stiffness and Mass Matrices . Program 43 : Equilibrium of Two-dimensional Rigid-jointed Frames Program 44 ; Equilibrium of Three-dimensional Rigid- jointed Frames. Program 43 : Equilibrium of Two-dimensional Pin-jointed Trusses... : Program 46 ; Equilibrium of Three-dimensional Pin-jointed Trusses : Program 47 : Plastic Analysis of Two-dimensional Frames Program 48 : Stability Analysis of Two-dimensional Frames Program 49 : Equilibrium of Rectangular Plates in Bending 104 107 112 47 121 124 127 134 138Chapter 5 Static Equilibrium of Linear Elastic Solids... . . $0 SA Introduction Program 50 Program $.1 Program $2 Program 53 Program 5.4: Program $5 : Program 5.6 Program 5.7: Program $8 Program 59 : Program $.10: References Plane Stress Analysis Using Three-node Triangles Plane Stress Analysis Using Six-node Triangles 20 Plane Strain Analysis Using Fifteen-node Triangles : Plane Strain Analysis Using Four-node | Quadrilaterals Plane Strain Analysis Using Eight-node Quadrilateral Plane Strain Analysis Using Nine-node Quadrilaterals . ‘ Axisymmetric Strain Analysis Using Four- | node Quadrilaterals . . ; Non-axisymmetric Strain of Axisymmetie Solids Using Eight-node Quadrilaterals . . : Three-dimensional Analysis Using Four-noded Tetrahedra ‘Three-dimensional Analysis Using Eight- noded Brick Elements . : ‘Three-dimensional Analysis Using Twenty noded Brick Elements Chapter 6 Material Nonlinearity . 2... ee we 60 61 62 63 64 65 66 67 Introduction Stress-Strain Behaviour | ol) ll. Stress Invariants, Failure Criteria Generation of Body-loads. | 1) le Visco-plasticity Initial Stress Corners on the Failute and Potential Surfaces Program 60: Program 6.1: Plane Strain Bearing Capacity Analysis— ‘Visco-plastic Method Using the Von Mises Criterion . ae Plane Strain Slope Stability Analysis Visco plastic Method Using the Moht-Coulomb Criterion prec 1st 142 143 143 144 150 154 157 161 165 169 172 178 182 186 191 192 192 195 197 19 200 201 202 203, 21Program 6.2 : Plane Strain Passive Earth Pressure Analysis—Initial Stress Method Using the Mohr-Coulomb Criterion Seec Program 6,3 : Axisymmetric ‘Undrained’ Analysis —Visco- plastic Method Using the Mohr-Coulomb Criterion Program 64 : Three-dimensional Elasto-plastic Analysis— Visco-plastic Method Using the Mohr Coulomb Criterion... 68 References... Gonos 5555568 Chapter7 Steady State Flow . Program 7.0 : Solution of Laplace's Equation Over a Plane Area Using Four-node Quadrilaterals. Program 7.1 : Analysis of Planar Free Surface Flow Using Four-node Quadrilaterals te 71. Reference Chapter 8 Transient Problems: First Order (Uncoupled) 8.0 Introduction we Seb oodg Program 80 : Solution of the Conduction Equation over a Rectangular Area . Program 8.1 : Solution of the Conduction Equation over a Cylinder Program 82 : Solution of the Conduction Equation over a Rectangular Area Using an Explicit Method Program 8,3 : Solution of the Conduction Equation over a Rectangle using an EBE Product Algorithm 8.1 Comparison of Programs 80, 8.2 and 8.3, : Program 84 : Diffusion~Convection Problem over a Rectangle—Transformed Analysis Program 85 : Diffusion-Convection Problem over a Rectangle—Untransformed Analysis 82. References . a Chapter 9 ‘Coupled’ Problems 2... 2. 9.0 Introduction... Program 90 : Solution of the Steady-state Navier-Stokes Equations over a Rectangle... . 9.1 Analysis of Soil Consolidation using Biot’s Theory . Program 9.1 : Biot Consolidation of a Rectangle in Plane Strain (Four-node Elements) Program 9.2 : Biot Consolidation of a Rectangle in Plane Strain (Eight-node/Four-node Elements) . 92 References... poudo oD ee - 18 224 - 232 237 239 239 245 250 . 250 - 250 256 259 = 262 + 26 . 267 2m 216 217 . 21 218 286 + 286 - 294 + 299Chapter 10 Eigenvalue Problems 100 tot Chapter 11 110 td Introduction Program 10\ Program 10.1 Program 10,2 : Program 103 Program 104 References Introduction, Program 11.0 : Program 11.1 : Program 11.2 : Program 11.3 : Program 11.4 ; Program 11.5 : References igenvalues and Eigenvectors of a Two- dimensional Frame of Beam-Column Elements—Consistent Mass Eigenvalues of a Rectangular Solid in Plane Strain—Four-node Quadrilaterals with Lumped Mass . Eigenvalues and Eigenvectors of a Rectangular Solid in Plane Strain—Eight- node Quadrilaterals with Lumped Mass . Eigenvalues and Eigenvectors of a Rectangular Solid in Plane Stress—Bight- node Quadrilaterals with Consistent Mass . : Eigenvalues and Eigenvectors of a Rectangular Solid in Plane Strain—Four- node Quadrilaterals with Consistent Mass Forced Vibrations. Forced Vibration of a Rectangular Solid in Plane Strain Using Eight-node Quadrilaterals, Lumped Mass, Modal Superposition Method Forced Vibration of a Rectangular Solid in Plane Strain Using Eight-node Quadrilaterals—Lumped or Consistent Mass, Implicit Integration by Theta Method Forced Vibration of a Rectangular Elastic Solid in Plane Strain Using Eight-node Quadrilaterals—Lumped or Consistent Mass, Implicit Integration by Wilson Theta Method . Forced Vibration of a Rectangular Solid in Plane Strain Using Eight-node Quadrilaterals—Lumped Mass, Complex Response Method. Forced Vibration of a Rectangular Elasto- plastic Solid in Plane Strain Using Eight-node Quadrilaterals—Lumped Mass, Explicit Integration Forced Vibration of an Elastic Solid in Plane Strain Using Four-node Quadrilaterals— Lamped or Consistent Mass, Mixed Explicit/Implicit Integration ", 300 304 308 312 35 319 320 320 320 325 330 333 337 345 + 350x Chapter 12 Analysis of Piles and Pile Groups... 2... 120. Introduction an Program 120,: t-2 Analysis of Axially Loaded Piles Using ‘Two-noded Line Elements. Program 121% p-y Analysis of Laterally Loaded Piles Using Two-noded Beam Elements Program 122 : Analysis of Vertically Loaded Pile Groups Using Two-noded Line Elements—t-z Spring forthe Soil and Mindlin's Equation for the Interactions Program 123 : Pile Drivability by Wave Equation—Two- node Line Elements—Implicit Integration in Time 121 References . Appendix 1 Consistent Nodal Loads Appendix 2 Plastic Stress-Strain Matrices and Plastic Potential Derivatives Appendix 3 Geometry Subroutines Appendix 4 Alphabetic List of Building Block Subroutines Appendix 5 Listings of Special Purpose Routines ‘Author Index Subject Index 351 1 381 . 38l . 389 - 365 372 378 379 382 386 401 415 466Preface to Second Edition Following the success of the First Edition, considerable modifications and improvements have been made. In retrospect it was probably a mistake to combine in one book a programming philosophy and a survey of other workers’ applications to geomechanics, so the latter has been dropped. However, the ‘geomechanics flavour remains, particularly in Chapters 6, 7, 8, 9 and 12 Although the prograraming philosophy remains intact, Chapter I reflects the rapid developments in the hardware and (to a lesser extent) software fields. Emphasis is given to what are generally known today as micro’ computers while the programming language has been updated from FORTRAN IV (FORTRAN 66) to FORTRAN V (FORTRAN 77), Chapter 2 has been extended to give a fuller treatment of structural applice ations and the equations of fuid flow including the Navier-Stokes equations. Inthe original Chapter 3, choice of storage strategy was rather limited, so more versatility has been introduced, in the form of a ‘profile’ (‘skyline’) option. Chapter 4 has been extended to embrace elastic, elasto-plastic and clastic stability analyses of ali commonly encountered skeletal structures, while in Chapter 5a much greater range of plane and solid finite elements is provided. In Chapter 6 more realistic examples of analyses of geotechnical problems — involving walls, foundations and slopes—have been incorporated. The original Chapter 7 dealt both with steady state and with transient flow problems. In the present edition, transient problems are removed to Chapter 8 while the steady state confined and unconfined flow problems analysed in Chapter 7 are more typical of those encountered in geotechnical engineering Practice. Chapter 8 also includes relatively recent developments of element-by-element algorithms for first order transient problems. ‘The original Chapter &, now Chapter 9, has been extended to deal with two typical coupled” problems—those arising from the (steady state) Navier-Stokes equations and those arising from Biot’s equations for transient analysis of porous elastic solids. ‘The original Chapter 9 has been subdivided, so that eigenvalue analysis (new Chapter 10) is divorced from the solution of second order transient problems (new Chapter 11), In the new Chapter 10 a solver based on the Lanczos method for large eigenvalue problems is demonstrated. In the new Chapter 11 second xixii order differential equations are integrated by a greater range of methods, for example mixed integration operators. ‘The final chapter (new Chapter 12) has been expanded to embrace the analysis of pile groups as well as of single piles. COMPUTER PROGRAMS All software (libraries and programs) described in this book are available from: Numerical Algorithms Group Ltd, NAG Central Office, Mayfield House, 256 Banbury Road, ‘Oxford OX2 7DE, UK Specify the medium (tape, disc) and format when applying. A small handling charge is made.CHAPTER 1 Preliminaries: Computer Strategies 10 INTRODUCTION Many textbooks exist which describe the principles of the finite clement method of analysis and the wide scope of its applications to the solution of practical engineering problems. Usually, little attention is devoted to the construction of thecomputer programs by which the numerical results are actually produced. It is presumed that readers have access to pre-written programs (perhaps to rather complicated ‘packages’) or can write their own. However, the gulf between understanding in principle what to do, and actually doing it, can still be large for those without years of experience in this field ‘The present book attempts to bridge this gulf. Its intention is to help readers assemble their own computer programs to solve particular engineering problems by using a ‘building block’ strategy specifically designed for computations via the finite clement technique. At the heart of what will be described is not a ‘program’ ‘ora set of programs but rather a collection (library) of procedures or subroutines which perform certain functions analogous to the standard functions (SIN, SQRT, ABS, ete) provided in permanent library form in all useful scientific ‘computer languages. Because of the matrix structure of finite element formula- tions, most of the building block routines are concerned with manipulation of matrices. ‘The aim of the present book is to teach the reader to write intelligible programs and to use them. Super-efficiency has not been paramount in program construction. However, all building block routines (numbering over 100) and all test programs (numbering over 50) have been verified on a wide range of computers. Their efficiency is also believed to be reasonable, given a tolerably Well written compiler. In order to justify the choice of high level language used (FORTRAN 77) and the way in whichitis asserted programs should be written, it is necessary briefly to review the development of computer ‘architecture’ and its relationship to program writing 1.1 COMPUTER HARDWARE In principle, any computing machine capable of compiling and running FORTRAN programs can execute the finite element analyses described in thisTable 1.1 ‘Approximate Wordlength Core storage cost 1985, Designation Example (its) (words) Virtual memory? LUK) Micro TBM PC 16 No 2k Super micro ICL PERQ 2 Yes 20K Mini vax 11/780 2 Yes 100K Mainframe CDC 7600 6 No 000K. Yes 10.000K Super mainframe CDC 10GF Cy book. In practice, hardware will range from ‘micro’ computers for small analyses and teaching purposes to ‘super’ computers for very large (especially non-linear) analyses. It is a powerful eature of the programming strategy proposed that the same software will run on all machine ranges. The special features of vector and array processors are described in a separate section (1 Table 1.1 shows that, in very general terms, five categories of useful hardware can be distinguished at the time of writing, graded using an ‘order of magnitude’ rule with respect to cost. While the trend in costs, certainly for minis’ and below, is downwards some such cost-based tabulation will always apply. At the lower end of the cost range, itis of course possible to code finite element analysis algorithms for very cheap machines, by employing machine code or an interpreter such as BASIC. The complete lack of portability of the former and the lack of a true subroutine structure inthe latter leads the authors to exclude such hardware from consideration. The user’s choice of hardware is a matter of accessibility and of cost. For example, in very rough terms, a “mini” computer (Table 1.1) has a processing speed about 30 times slower than that of a ‘mainframe’. Thus a job taking two minutes on the mainframe takes one hour on the mini, Which hardware is “better” clearly depends on individual circumstances. The main advice that can be tendered is against using hardware that is too weak for the task; that is the user is advised not to operate at the extremes of the hardware’s capability. 1.2. STORE MANAGEMENT. In the programs in this book it will be assumed that sufficient main random access core store is available for the storage of data and the execution of programs. However, the arrays processed in very large finite element calculations cean be of size, say, 10,000 by 200. Thus a computer would have to have a main store of 2 x 10° words to hold this information, and while some such computers exist, they are still comparatively rare. A much more typical core size is still of the order of 1 x 10° words. Two distinct strategies are used by computer manufacturers to get round the problem of limited random access core. In the first (simplest and earliest)3 technique different ‘Tevels' of store are provided with transfer of data to core from various secondary storage devices such as discs and tape being organized by the programmer. The CDC 7600 is a computer that used this strategy. Programs are executed in a small fast store while data can be stored in a medium- {ast ‘large core memory’. Further back-up storage is available on a big dise but of ‘course when the programmer decides to transfer information to and from dise the central processor would be idle, A job calling for disc transfers is therefore automatically swapped out of fast memory, and is replaced by the next job waiting in the queue. The second strategy removes store management from the uscr’s control and vests it in the system hardware and software. The programmer sees only a single level of store of very large capacity and his information is moved from backing, store to core and out again by the supervisor or executive program which schedules the flow of work through the machine. This concept, namely of a very large virtual core, was frst introduced onthe ICL ATLAS in 1961. Itcan be seen from Table |.1 that this strategy of store management is now very common. The ATLAS and most subsequent machines of this type also permit multi- programming in that there can be several partly completed programs in the core at any one time. The aim in any such system is to keep the hardware (CPU, peripherals, etc) fully utilised. Clearly it is necessary for the system to be able to translate the virtual address of variables into a reat address in core. This translation usually involves a complicated bit-pattern matching called ‘paging’. The virtual store is split into segments or pages of fixed or variable size referenced by page tables, and the supervisor program tries to ‘learn’ from the way in which the user accesses data in order to manage the store in a predictive way. However, store ‘management can never be totally removed from the user's control. It must always bbe assumed that the programmer is acting in a reasonably logical manner, accessing array elements in sequence (by rows or columns as organised by the compiler), I'the user accesses a virtual store of 10° words ina random fashion the paging requests will readily ensure that very little execution of the program can take place. 13 VECTOR AND ARRAY PROCESSORS Tnlater chapters it will be shown that finite clement computations consis largely ‘of handling matrices and vectors—multiplying them together, adding, transpos- ing, inverting and so on. It has seemed appropriate to computer manufacturers, therefore, to provide special hardware capable of doing vector or matrix operations at high speed. This hardware can be in the form of central processing, vector registers such as on the Cray and CDC Cyber 205 or in peripheral form such as the ICL 2900 DAP concept. Other array processors can also be connected as peripherals, one such being the Floating Point Systems FPS264, and this is a very likely future development, The drawback associated with specialised hardware of this kind is that4 specialised programs usually have to be written to exploit it, This detracts from such important software considerations as program portability and machine independence. However, at the very least, software can be written to pre-process standard high level language code to sec if it can be re-organised to make best use of a particular machine's vector capabilities. Future trends in this direction are difficult to anticipate, It may be that, ultimately, array processing will become so fast and cheap, even for very large arrays, that much of the special purpose software which has been written to date, for example to process arrays with various sparsity patterns, will become redundant. However, such a trend would be entirely consistent with the philosophy of the present book, and would merely necessitate a change of emphasis, At present, though, only standard, completely portable high level languages will be considered in the following sections on software. 14 SOFTWARE, Since all computers have different hardware (instruction formats, vector capability, etc) and different store management strategies, programs that would make the most effective use of these varying facilities would of course differ in structure from machine to machine, However, for excellent reasons of program portability and programmer training, engineering computations on all machines are usually programmed in “high level’ languages which are intended to be ‘machine-independent. The high level language (FORTRAN, ALGOL, ete) is translated into the machine order code by a program called a ‘compiler. Tt is convenient to group the major high level languages into those which are FORTRAN-like (for example PASCAL) and those which are ALGOL-like (for example ADA), The major difference between these groups i that the latter allow dynamic management of storage at run-time whereas the former do not. This difference has immediate implications for store management strategy because clear that unless a computer permits genuine multi-programming and aone-level store in which programs are dynamically expanding and contracting there is not much point in having a language that allows these facilities. Thus, on CDC 7600 machines, FORTRAN was an obvious choice of language whereas on other computers there is greater room for manoeuvre. Dynamic store management seems to be a uniform trend for the future, so one would expect ALGOL-like languages to become more and more attractive. ‘There are many texts describing and comparing computer languages, often from the point of view of computer scientists, ¢g. Barron (1977). What matters in the context of the present book is not syntax, which is really remarkably constant among the various languages, but program structure. The following example illustrates the two crucial differences between ALGOL-ike and FORTRAN-like languages as far as the engineer-programmer is concerned. ‘Suppose we want to write a program to multiply two matrices together, a task that occurs hundreds of times during execution of a typical finite clementprogram. In ALGOL 60 the program might read as follows: Degin integer J, m, n statements which read or compute I,m, n: Degin array a{I:1,1:mi, B[1-m,1:n], e[1:, 1:0} integer i,j,k; real x; statements which initialise a,b; for i:=1 step } until ! do for j:=1 step J until n do begin x for k:=1 step 1 until m do xm xt afisk] bk end; eli,fls=x statements which do something with ¢ end end while in FORTRAN 77 it might be written: INTEGER I, J, K, L, M, N REAL X REAL A(20, 20), B(20,20), C(20,20) Statements which read or compute L, M, N Statements which initialise A, B DO! 0 DO2K=1,M 2 X=X+A(LK)«BK,)) CU) =X 1 CONTINUE Statements which do something with C END Note that the declarations of the simple integers and the real in lines 1 and 2 of, this program are unnecessary. ‘Comparison of the two codes will casily show that the syntaxes of the two languages are very similar, eg. compare the ALGOL for i:=1 step 1 until | do with the corresponding FORTRAN DO1T=1,L or the assignment statements xrax+ali,k]*b Uk:6 and X=X+AGK)eBKD, ‘Therefore, anyone familiar with one syntax can learn the other in a very short time. However, consider the program structures. ‘The main difference is that in ALGOL, the array dimensions J, m,n do not need to be known at compile time whereas in FORTRAN they do. Indeed, our FORTRAN program is uscless if we ‘want to multiply together arrays with dimension greater than 20. Iti this lack of ‘dynamic’ storage which is FORTRAN’s major limitation, It can partially be alleviated by various devices such as allocating all storage to a large vector and keeping track ofthe starting addresses ofthe various arrays held inthe vector, but one can never get away without a dimensioning statement of some kind within a FORTRAN program. “This has various bad effects, among which isa tendency for FORTRAN users to declare their arrays to have the maximum size that will ever be encountered, subject to the limitation ofthe core size of the machine, Since most large machines allow multi-programming this can be inefficient. 15 PORTABILITY OF ‘SUBROUTINES’ OR ‘PROCEDURES’ ‘The second major deficiency of FORTRAN becomes apparent when we wish to construct a program to multiply together matrices of varying size within the same program, This sort of activity occurs very frequently in finite element programs. ‘To avoid the repetition of sections of code within a program, most languages allow these sections to be hived off into self-contained sub-programs which are then called” by the main program. These are termed ‘procedures’ in ALGOL and ‘subroutines’ in FORTRAN and their use is central to efficient and readable programming, For example, it would be very tedious to include the code required to compute the sine of an angle every time it is required, so ALGOL and FORTRAN provide a special, completely portable routine (actually a function) called SIN, This is held in a permanent library and is accessible to any user program. Why not do the same thing with matrix multiplication? Tn ALGOL the technique is straightforward and encouraged by the structure of the language. For example, an ALGOL 60 matrix multiplication procedure might be as follows: procedure matmult (a,,¢, m,n); value I, m,n; array ab, ¢; integer I, m, n; begin integer i,j, k: real x; step ! until n do begin x:=.0 for k:= 1 step I until m do xrext afk] +bLb fl:end oli, ls end Suppose we now wish to multiply the (! x m) array a by the (m x n) array b to give the (I n) array c as before, and then to multiply c by the (n xl) array d to yield the (1 I) array e. The ALGOL 60 program would be as follows: begin integer |, m,n; procedure matmult (a,b, ¢,1.m,nly statements which read or compute I, m, n; begin array af1:l,1:ml, BEI:m, 1:7), eCUeh tem), don 0, eC: statements which initialise a, by rmatmult (a,b,c 1,m4n) rmatmult (c,d.e, statements which do something with ¢ end end Inthe above, the procedure matmult is assumed to be available in a user library in exacily the same way as LOG, COS, etc, are available in the permanent library. The program structure can be scen to be simple and logical, reflecting exactly what the user wants to do and requiring a minimum of information about array dimensions to do it. Let us now try to do the same thing in FORTRAN 77. Reasoning ftom out previous program we write: SUBROUTINE MATMUL (A, B,C,L,M,N) REAL A(20, 20), B(20, 20}, C(20,20), X INTEGER I, J, K, L, M, N, DOI x=0 DO2K=1,M 2 X=X+A(LK)eBKD, CUJ)=X 1 CONTINUE RETURN. END. ‘Now the dimensioning statement in line 2 of this subroutine is merely a device for passing the addresses in memory of the arrays locally called A, B, C back and forth between the subroutine and the main program. Since the FORTRAN subroutine only needs to know a starting address and a column size, we could as well have written REAL A(20, 4}, B(20, 6), C(20,»), XMore importantly, since in FORTRAN the local A, B, C are mapped onto the arrays of the main program, it can be seen that our subroutine is useless, unless arrays processed by it have a first dimension of 20. This is an impossible restriction akin to saying we can only work out the sines of angles less than 20° using SIN. ‘A way round this deficiency in FORTRAN is to make the subroutine completely portable at the expense of including extra parameters in the call, associated with each array. These parameters, which are constants, contain on entry to the subroutine the column sizes of the particular arrays to be multiplied at that time as declared in the main program. Thus our portable FORTRAN subroutine becomes SUBROUTINE MATMUL (A, IA, B, IB, C, IC, L, M, N) REAL A(IA, ¢), B(IB,*), C(IC #), X INTEGER IA, IB, IC, L, M,N, I,J, K 2 X=X+A(,K}«B(KD cuy=X 1 CONTINUE RETURN, END We could now construct our general matrix multiplication program asfollows, assuming MATMUL to be available in a user library: INTEGER L, M, N, IA, TB, IC, 1D, IE REAL A(20, 5}, B(S, 15}, C(20, 15), D(15,20), E(20,20) DATA 1A/20), 1B/5), 1C)/20/, 1D/15), 1B,/20/ Statements which read or compute L, M, N Statements which initials A, B CALL MATMUL (A,1A,B,IB,C,IC, L, M,N) : : CALL MATMUL (G.IC,D,1D,E,1E,L,N,L) Statements which do something with E END We have now succeeded in making FORTRAN look reasonably like ALGOL. ‘The differences lie only in the need to dimension arrays in the main program and to pass the extra column size parameters in the subroutine calls. In fact we will only be concerned in our finite element programs with adjusting the sizes ofa few large arrays at run-time, and so the inconvenience can really be restricted to changing a line at the beginning of the main program. “This is most easily done by the PARAMETER statement in FORTRAN 77. ‘Using this facility the second and third lines of the preceding program can be replaced by:PARAMETER (1A = 20, JA =5, IB=5, JB=15, IC = 20, JC=15, *1D = 15, JD =20, IE =20, JE =20) REAL A(IA,JA), B(IB.JB}, CIC, JC), DAD, JD), E(IE, JE) so that only the PARAMETER line need ever be changed. This used in this book the method 1.6 LIBRARIES ‘The concept that will be described in subsequent chapters is one of having a user library containing completely portable subroutines such as MATMUL, which was described in the previous section. The idea is not new, of course, and many such libraries already exist, notably the NAG mathematical subroutine library (Ford and Sayers, 1976). The present library is completely compatible with the NAG mathematical subroutine library and uses the same parameter passing convention, but it differs from the NAG mathematical subroutine library in two ways. Firs, it contains only routines used in finite element analysis, and second it uuses mnemonic names such as MATMUL to describe its operations (an equivalent routine in the NAG mathematical subroutine library is called FOICKF). This greatly enhances program readability. The library routines which are the “building blocks’ from which the finite clement programs are constructed are described in Chapter 3. Some of them may bbe considered as ‘black box’ routines, that is the user need not understand or see them (how many users know how SIN(X)is computed”). Some of these routines are not even listed in the book but the source code is available on tape or floppy disc in all the usual formats. Other routines are, however, specifically related to the descriptions of finite element operations, and these routines are listed for reference purposes in Appendix 5, and of course on tape also, ‘The test programs which are built up from the library routines are described in Chapters 4 to 12 and also listed there. Interested readers can substitute their own routines at any point, because of the modular nature of the whole system. Such routines could be in machine code for greater efficiency or could be linked to the use of an array processor. For example, the Floating Point Systems Array Processor can look to the user like a set of library subroutine calls which could be substituted for the FORTRAN subroutines if required. ‘The programs described in this book arc all assumed to be processed ‘in core’. "Naturally, one could construct out of core versions ofthe library routines, which Process large arrays for usc on non-virtual machines. 1.7, STRUCTURED PROGRAMMING The finite element programs which will be described are strongly ‘structured? in the sense of Dijkstra (1976), The main feature exhibited by our programs will be seen to be a nested structure and we will use representations called ‘structure charts’ (Lindsey, 1977) rather than flow charts to describe their actions.10 ‘The main features of these charts are: (@ The Block DOTHS DOTHAT DOTHEOTHER ‘This will be used for the outermost level of each structure chart. Within a block, the indicated actions are to be performed sequentially. (ii) The Choice Anawert | Answer? | Answer3 ‘ACTION’ | ACTION2 | ACTIONS ‘This corresponds to the if...then...else if...then...end if kind of cohstruct, ii) The Loop ‘This comes in various forms, but we shall usually be concerned with the‘DO' loop ACTION TOBE REPEATED TIMES In particular, the structure chart notation encourages the use only of those GOTO statements that are safe in a structured program. Using the notation, our matrix multiplication program would be represented as follows:ul arraysaand row x column products faxb) Do something withe The nested nature of a typical program can be seen quite clearly. 18 OTHER LANGUAGES ‘New programming languages appear from time to time, usually at the instigation ofcomputer scientists rather than engineers. For example, PASCAL is popular as 8 teaching language and has a simple structure enabling compact compilers to be written. In PASCAL 3 our matrix multiply algorithm becomes: TYPE ARR = ARRAY [INTEGER, INTEGER] OF REAL; PROCEDURE MATMUL (VAR A, B, C: DYNAMIC ARR; L, M, N: INTEGER) VAR I, J, K: INTEGER; X:REAL; BEGIN FOR I:=1 TOL DO FOR J:=1 TO N DO BEGIN X:= 00;12 FOR K:=1 TO M DO XX +ALL K]+BIK, J], CUI Xs END; END; VAR L, M, N: INTEGER; ‘A:(L-10, 1.10] OF REAL; B:[1.10, 1.10] OF REAL; C:[1..10, 1.10] OF REAL; BEGIN Statements which read L, M,N, A, B, C MATMUL (A,B,C,L,M,N); Statements which do something with C END. ‘Thus the language has some ALGOL-like features in that, within the procedure, arrays A, B and C are quasi-dynamic and do not need to be dimensioned. However, there is no real concept of a stack pointer and in the main program A, B and Creed to be assigned fixed dimensions. The PASCAL language does possess. the builtin procedures (functions) HIGH and LOW which enable bounds of arrays to be inspected: e.g. HIGH(A, 1) would deliver the upper bound of the first dimension of A. However, because of the lack of a truly dynamic structure, these cannot be used to simplify the procedure call statements. ‘A language in which this is possible is ALGOL 68. It has an extremely flexible structure and points the way towards the probable shape of future languages. Our matrix multiply program becomes: (proe maimul = ([,] real a, b, ref [,] real ¢) void: (int [= 1 upb a, m= upb b, n=2 upb; lo int i,j, ki for i to Ido for j to n do loc real x for k to m do xt malik] +bUk fled: efi, J]:=x od od ¥ . foe int 1, m,n; read ((I,m,n)); Joe [1:i, sm} real a, loc [1:m, 1:1] real b, loe (1:1, 1:7] real ¢: read ((a,b)); ‘matmul (a,6,¢); print (¢) ) ‘The procedure calling parameters have been reduced to three (compared with FORTRAN's nine) and the language has many more useful features. For3 example, operations between matrices can be defined so that procedure calls can be dispensed with and operation c= ax b would implicitly mean matrix multiply. One can therefore anticipate languages specially adapted for matrix manipulation. 19 CONCLUSIONS Computers on which finite element computations can be done vary widely in their size and architecture. Because of its entrenched position FORTRAN is the language in which computer programs for engineering applications had best be written in order to assure maximum readership and portability. (There are inevitably minor differences between machines, for example in their handling of free format input which is used exclusively in the test programs, but these are not overwhelming.) However, FORTRAN possesses fundamental limitations which can only partially be overcome by writing a structure of FORTRAN which is as, close to ALGOL as can be achieved. This is particularly vital in subroutine calling, Using this structure of FORTRAN, a library of subroutines can be created which is then held in compiled form on backing store and accessed by programs in just the way that a manufacturer's permanent library is. Using this philosophy, a library of over 100 subroutines has been assembled, together with some 50 example programs which access the library. These programs are written in a reasonably ‘structured’ style, and tapes or discs of library and programs will be supplied to readers on request. Versions are at present available for all the common machine ranges. The structure of the remainder ofthe book is as follows. Chapter 2 shows how the differential equations governing the behaviour of solids and fluids are programs are written in a reasonably ‘structured! style, and tapes or discs of library and programs will be supplied to readers on request. Versions are at present available for all the common machine ranges. Chapter 3 describes the sub-program library and the basic techniques by which main programs are constructed to solve the equations listed in Chapter 2 ‘The remaining Chapters 4 to 12 are concerned with applications, partly in the authors’ field of geomechanies. However, the methods and programs described are equally applicable in many other fields such as structural mechanics, water resources, electromagnetics and so on. Chapter 4 leads off with static analysis of skeletal structures and plate problems. Chapter $ deals with static analysis of linear solids, whilst Chapter 6 discusses extensions to deal with material non- linearity. Chapter 7 is concerned with problems of fluid flow in the steady state while transient states with inclusion of transport phenomena (diffusion with advection) are treated in Chapter 8. In Chapter 9, coupling between solid and. fluid phases is treated, with applications to consolidation processes. A second type of ‘coupling’ involves the Navier-Stokes equations. Chapter 10 contains programs or the solution of steady state vibration problems, involving the deter- mination of natural modes, by various methods. Integration of the equations4 ‘of motion in time by direct integration, by modal superposition and by ‘complex response’ is described in Chapter 11. Finally, Chapter 12 deals with problems of pile capacity and driveability, for single piles and including group effects. In every applications chapter, test programs are listed and described, together with specimen input and output 1.40 REFERENCES Barron, D. W.(1977) An Introduction tothe Study of Programming Languages. Cambridge University Press. Dijkstra, E. W. (1976) A Discipline of Programming. Prentice-Hall, Englewood Cliff, New Tersey, Ford, Band Sayers, D. K. (1976) Developing a single numerical algorithms library for diferent machine ranges. ACM Trans. Math. Software, 2, 115. Lindsey, C. H. (1977) Structure charts a structured alternative to flow charts. SIGPLAN Notices, 12, No. 11, 36CHAPTER 2 Spatial Discretisation by Finite Elements 2,0 INTRODUCTION ‘The finite element method is a technique for solving partial differential equations by first discretising these equations in their space dimensions. The discretisation is carried out locally over small regions of simple but arbitrary shape (the finite elements). This results in matrix equations relating the input at specified points in the elements (the nodes) to the output at these same points. In order to solve equations over large regions, the matrix equations for the smaller sub-regions are usually summed node by node, resulting in global matrix equations. The method is already described in many texts, for example Zienkiewicz (1977), Strang and Fix (1973), Cook (1981), Connor and Brebbia (1976) and Rao (1982), but the principles will briefly be described in this chapter in order to establish a notation and to set the scene for the later descriptions of programming techniques. 21 ROD ELEMENT Figure 2.1(a) shows the simplest solid element, namely an clastic rod, with end nodes 1 and 2. The clement has length L while u denotes the longitudinal displacements of points on the rod which is subjected to axial loading only If Pis the axial force in the rod at a particular section and F is an applied body force with units of force per unit length then P cA= Edm BAM ey ‘and for equilibrium from Figure 2.1(b), op ethno @2) Hence the differential equation to be solved is EATS +F=0 23) oF 1516 displacement c- / a _/dody force F (a) ee od + 7 +L —I oer e — —e Po dh x wo) ps 38 = rAd dtu LA ci ( —_— =p. 228 sae Figure 2.1 Equilibrium of a rod element (Although this isan ordinary differential equation, future equations, for example (2.13), preserve the same notation.) In the finite element technique, the continuous variable w is approximated in terms ofits nodal values, 1, and u,, through simple functions ofthe space variable called ‘shape functions’. That is waeNyuy + Nats or usiM, waft} 4) {In simple examples the approximate equality in (24) could be made exact by vetting x M= (: = : Qs) ifthe true variation of wis linear. In general N, and N, will be of higher order or ‘more linear element subdivisions will be necessary. When (24) is substituted in (23) we have BASIN, waft} rao e6 so that the partial differential equation has been replaced by an equation in the Giscretised space variables w, and u,. The problem now reduces to one of finding, ‘good’ values for u, and tl, ‘Many methods could be used to achieve this, and Crandall (1956) discussesW collocation, subdomain, Galerkin and least squares techniques. Of these, Galerkin’s method, e.g. Finlayson (1972) is the most widely used infinite element ‘work. The method consists in multiplying equation (26) by the shape funetions in turn and integrating the resulting residual over the element, Thus LON) a 7 rey, Nea Son, Nal} Vrdeed 27 [jordin safsjan irae on ‘Note that in the present example, double differentiation ofthe shape functions would cause them to vanish and yet we know the correct shape may not be of higher order than linear. This difficulty is resolved by applying Green's theorem (integration by parts) to yield typically [ritioce [2 Mace nomday tem ir aioe 25 Hence (2.7) becomes: Ox ox ON, ON; ox ox | LON aN, ON: aft} {nef 9 @9) ed On evaluation of the integrals, pe =| 2.10) ‘The above case is for a uniformly distributed force F acting along the rod. For the case in Figure 2.1(a) where the loading is applied only at the nodes we have af HSH on 1 which are the familiar ‘stiffness’ equations for a uniform prismatic rod. In matrix notation we may write KMu=F 212) Where KM is called the ‘element stiffness matrix’. 2.2 ROD INERTIA MATRIX Consider now the case of an unrestrained rod in longitudinal motion, Figure 2.1(¢ shows the equilibrium ofa segment in which the body force is now18 siven by Newton's law as mass times acceleration. Ifthe mass per unit volume is », the differential equation becomes eu ae Bana OA ae 213) (On discretising u by finite elements as before, the first term in (2.13) clearly leads again to KM. The second term takes the form TIN) uy -{ {Rhea Sem, wal {it}ex ey HENAN, NN2] 4. & fur ~o[ [ax mate ors Evaluation of integrals yields $4] 2 fu oat} jhe fe} - or in matrix notaion -wm at oF where MM is the ‘clement mass matrix’ or ‘element inertia matrix’. Thus the full matrix statement of equation (2.13) is ae KMu+MM<> at @17) ‘Note that MM as formed in this manner is the ‘consistent’ mass matrix and differs fror: the ‘lumped’ equivalent which would lead to 4pAL terms on the diagonal with zeros off-diagonal. 23. THE EIGENVALUE EQUATION Equation (2.17) is sometimes integrated directly (Chapter 11) but is also the starting point for derivation of the eigenvalues of an element or mesh of elements. ‘Suppose the elastic rod element is undergoing free harmonic motion. Then all ‘nodal displacements will be harmonic, of the form usasin(ot +) (218) where a are amplitudes of the motion, « its frequencies and y its phase shifts. ‘When (2.18) is substituted in (2.17) the equation KMa—o?MMa=0 219) is obtained, which can easily be rearranged as a standard eigenvalue equation. Chapter 10 describes solution of equations of this type.19 24 BEAM ELEMENT ‘As a second one-dimensional solid clement, consider the slender beam Figure 2.2. The end nodes 1 and 2 are subjected to shear forces and moments which result in translations and rotations. Each node, therefore, has two degrees of freedom. ‘The element shown in Figure 22 has length L, flexural rigidity EJ and carries a uniform transverse load of q pet unit length. The well known equilibrium equation for this system is given by a (2.20) Again the continuous variable, win this case, is approximated in terms of discrete nodal values, but we introduce the idea that not only w itself but also its derivatives can be used in the approximation. Thus we choose to write wetNy Na Ns Ns] 221) where 8, = dw/8x at node 1 and so on. In this case equation (2.21 can often be made exact by choosing the cubic shape functions: Nya p(t? 3s +204) Naps = 21x? +3) 1 (2.22) Ny=fpOLet—20) Nea fet-be) Ke A mC owe Sy Fares -—t — (Psa, oisiacenents Figure 22. Slender beam element20 Note that the shape functions have the property that they or their derivatives in this case equal one at a specific node and zero at all others. Substitution in (2.20) and application of Galerkin’s method leads to the four clement equations: % wy Ny f Nah ey IN, NA Ny Nad” of 3 dx (223) Jo) vy Pan EN ie Ne 6 Ng Again Green’s theorem is used to avoid differentiating four times; for example ay, ON, ON, aN, @N, fu Sdte=- ian = [FA Mae + nels tems (2.24) Hence assuming that E/ and q are not functions of x, (2.23) become ws ™, af SMM get jm,2,34 4 =f Nh dx @25 O Ng Evaluation of the integrals gives zw 6 2 6) L B PR PG 2 4.6 21/1], B ql L 2 L 1 ary EL B eae (2.268) Symmetrical 12-6] | w, 2 v7 B -Bl lm ‘ 4 B a n which recovers the standard slope-deflection equation for beam elements, ‘Theabove case is fora uniformly distributed load applied to the beam. For the case where loading is applied only at the nodes we have 12 6 2 6 Bp ogee #\{™) [* + -§ z a] |My a oo - 2260) re eet (aca ead a [Symmetrical tj led lw which represents the element stiffness relationship.2 Hence, in matrix notation we have KMw=F @2n Beam-column elements, in which axial and bending effects are superposed from (2.26b) and (2.11), are described further in Chapter 4. 25 BEAM INERTIA MATRIX If the element in Figure 2.2 were vibrating transversely it would be subjected to an additional restoring force — (0 w/ét?). The inertia or mass matrix, by analogy with (2.15), is just NUN NiNa NNy NANG ™ “1 NaN, NaNa NaNs NaNaly. 2 J 0, ~oal) NN, NGNy NGN, NGNG|* GP] w,p 278) NeN, NaNa NaNy NaN, a: Evaluation of the integrals yields 156 22L $4 = 130] py pat] a se 32 | o Jo, 229 420 156 -22L| oF Jw, (2.29) Symmetrical 4L* 6; Inthis instance, the neglect ofthe consistent mass terms leads to large errors in the prediction of beam frequencies as shown by Leckie and Lindburg (1963) Superposition of (2.29) and (2.16) is required in the dynamic analysis of framed structures and this is described further in Chapter 10. 26 BEAM WITH AN AXIAL FORCE Ifthe beam element in Figure 22 is subjected to an additional axial force P (Figure 2.3), a simple modification to (2.20) results in the differential equation aw, pow ElzatP aad (2.30) where the positive sign corresponds to a compressive axial load and vice versa. ro 1 b——. —4 Figure 23 Beam with axial force2 Finite element discretisation and application of Galerkin's method leads toan. additional matrix associated with the axial force contribution “ON ON, a Se aes 231) Evaluation of these integrals yields for compressive P 36 36 Toy -EZ 3] (w P 4c -3 -L| | +30 % | |, (232) fe 2 Symmetrical aL] | o, If this matrix is designated by KP the equilibrium equation becomes (KM — KP)w =F 233) Buckling of a member can be investigated by solving the eigenvalue problem when F = 0, by increasing the compressive force on the element (corresponding to KP in 2.33) until large deformations result or in simple cases by determinant search, Program 12.1 uses the KP matrix to assess the effect of axial loadingon the response of laterally loaded piles. Equations (2.32) and (2.33) represent an approximation of the approach to modifying the element stiffness involving stability functions (Lundquist and Kroll, 1944; Horne and Merchant, 1965). The accuracy of the approximation depends on the value of P/P, for each member, where Pr is the Euler load. Over the range — 1 < P/P, <1 the approximation introduces errors no greater than 7% (Livesley, 1975) For larger positive values of P/Pp, however, equation (2.33) can become highly inaccurate unless more element subdivisions are used. Program 4.8 therefore uses the stability function approach to modify the element stiffness matrices in analysing the stability of plane frames. 2.7 BEAM ON ELASTIC FOUNDATION In Figure 24 a continuous elastic support has been placed beneath the basic clement. If this support has stiffness k (force/length?) then clearly the transverse load is resisted by an extra force + kw leading to the differential equation er thw mg 234) By comparison with the inertial restoring force — p.A(@*w/ét®)it will be apparent that application of the Galerkin process to (2.34) will result in a foundationB e A UTititiiiat aS, Foundation stiffness k Figure 24 Beam on continuous elastic foundation stiffness matrix that is identical to the consistent mass matrix apart from the multiple + kL in place of —pAL. An example of a consistent finite element solution to (2.34) is given by Program 4.2 in Chapter 4. A lumped mass’ approach to this problem is also possible by simply adding the appropriate spring stiffness to the diagonal terms of the beam stiffness matrix and Chapter 12 describes programs of this type. 28 GENERAL REMARKS ON THE DISCRETISATION PROCESS Enough examples have now been described for a general pattern to emerge of how terms in a differential equation appear in matrix form after discretisation. Table 2.1 gives a summary, N being the shape functions. In fact, first order terms such as u/@x have not yet arisen. They are unique in Table 2.1 in leading to matrix equations which are not symmetrical, as indeed would be the case for any odd order of derivative. We shall return to terms of this type in due course, in relation to advection in fluid flow. Table 21 Semi-discretisation of partial differential ‘equations Term appearing Typial orm ia in difectia' te amen santos sate eguanon " fora tu oN, ie fries ou 2M, ON oe Ox Ox | ca ae oe ou fe an,u 29 ALTERNATIVE DERIVATION OF ELEMENT STIFFNESS Instead of working from the governing differential equation, element properties can be derived by an alternative method based on a consideration of energy. For example, the strain energy stored due to bending of a very small length of the beam element in Figure 2.2 is (235) where M is the ‘bending moment’ and by conservation of energy this must be equal to the work done by the external loads g; thus 4W=}qwox (236) To proceed, we discretise the displacements in terms of the element shape functions (2.21 and 2.22) to give w= Nyy + N38, +Nywy NO, = Nw 237) ‘The bending moment M is related to w through the ‘moment-curvature’ expression aw M=-EIE> or M=DAw (2.38) ‘where D is the material property EI and A is the operator ~ 6/63, Writing (235) in the form au -}(-25) an re) we have 6U =HAwy™M dx 240) Using (2.37) and (2.38) this becomes 6U = {ANw)'DANwx w™(AN)" DANw8x eat) The total strain energy of the element is thus u =f wT(ANTDANw dx (242) The matrix AN is usually written as B, and since w are nodal values and therefore constants er 243) o25 ‘Similar operations on (2.36) lead to the total external work done and hence the stored potential energy of the beam is given by n=u-w <0" [nar wg fines (244) A state of stable equilibrium is achieved when x is a minimum with respect to all w. That is am _ tee 7 guts |, BDBdew—a | Nidx=0 (2.45) [fmronaswaa nts 246) which is simply another way of writing (2.25). ‘Thus we see from (2.27) that the element stiffness matrix KM can be written in the form KM= f B'DBdx (247) which will prove to be a useful general matrix form for expressing stiffnesses ofall solid elements. The computer programs for analysis of solids developed in the next chapter use this notation and method of stiffness formation. 210 TWO-DIMENSIONAL ELEMENTS: PLANE STRESS AND STRAIN ‘The elements so far described have not been true finite elements because they have been used to solve differential equations in one space variable only. Thus the real problem involving two or three space variables has been replaced by a hypothetical, equivalent one-dimensional problem before solution. The elements ‘we have considered can be joined together at points (the nodes) and complete continuity (compatibility) and equilibrium achieved. In this way we can often ‘obtain exact solutions to our hypothetical problems and these solutions will be unaffected by the number of elements chosen to represent uniform line segments This situation changes radically when problems in two or three space dimensions are analysed. For example, consider the plane shear wall with ‘openings shown in Figure 2.5(a). The wall has been subdivided into rectangular elements of side a x b of which Figure 2.5(b) is typical. These elements have four corner nodes 50 that when the idealised wall is assembled, the elements will only be attached at these points. the wall can be considered to be of unit thickness and ina state of plane stress, see Timoshenko and Goodier (1951), the equations to be solved are the following:6 t Figure 2.5. (a) Shear wall with openings. (b) Typical rectangular element with four comer nodes @ Equilibrium 00, x By tty , 26, 48) Fgh where o,,0, and r,, are the only non-zero stress components and F,., F, are body forces, per unit volume, (i) Constitutive (plane stress) (Generalised Hooke) oh (as) where E is Young's modulus, v is Poisson's ratio and ¢,,e, and y,, are the independent small strain components. i) Strain-displacement a «) |x 0 a| fu l-fo aie) a fi fee Yo) Lay ox, where u,v are the components of displacement in the x,y directionsPu Using the notation of the previous section with A as the strain-displacement operator and D as the stress-strain matrix, these three equations become Max -F o=De (2.51) ende where a ox tv 0 efi} anjo yt oT esx a 0 ay ox We shall only be concerned in this book with ‘displacement’ formulations in which @ and ¢ are eliminated from (2.51) as follows: Mo=-F A'De= -F (2.53) A'DAe=-F Writing out (2.53) in full we have @u TD rey @ 2 a tay? Iavou t+y Hp which is a pair of simultaneous partial differential equations in the continuous space variables u and v. ‘As usual these can be solved by discretising over each element using shape functions w=[N; Nz Ns Ng) a =No (2.55) and vm IN, Nz Ny Neli'?}=Nv 2.56) where in the case of the rectangular element shown in Figure 2.5(b) the N,8 functions were first derived by Taig (1961) to be (257) ‘These result in linear variations in strain across the element which is sometimes called the ‘linear strain rectangle. Discretisation and application of Galerkin’s method (Szabo and Lee, 1969), using Table 2.1, leads to the stiffness equations for a typical element: (@N, aN), 1-vaN,Nj\/ ON, aN), 1- ON, aN, ) E ff \ ax ox TD By oy)\" ax dy 2 ey =) Jo do] (2M IML =n eNe A) (ANN, Ln aNLENY \"ay ax * 2 ox dy ay ty * 2 dx ox a) _ fe sdydx, j= 234 fh efit (2.58) g kMe=F Evaluation of the first term in the stiffness matrix yields E(b 1-va KM. 7a ($4575) om and so on In the course of this evaluation, integration by parts now involves integrals of the type i an ldxdy = ~ [ [Be Reavay + hu iras (2.60) ‘where lis the direction cosine ofthe normal to boundary S and we assume thatthe ‘contour integral in (2.60) is zero between elements. This assumption is generally, reasonable but extra care is needed at mesh boundaries. Only if the elements become vanishingly small can our solution be the correct one (an infinite number of elements) Physically, in a displacement method, it is usual to satisfy compatibility everywhere in a mesh but to satisfy equilibrium only at the nodes. It is also possible to violate compatibility, but none of the elements described in this book does.2 2.11 ENERGY APPROACH ‘As was done in the case of the beam clement, the principle of minimum potential energy can be used to provide an alternative derivation of (2.58). The element strain energy per unit thickness is Us [fieteaxay ait [faxrpanae dyr =r | [prpaxayr @61) where A and D are defined in (2.52) and N,N; Ny Ne 0 0 0 0 Nol /eiccsonioi wim Nims (262) Thus we have again for this element kM= fm DBdxdy 263) which is the form in which it will be computed in Chapter 3. Exactly the same expression holds in the case of plane strain, but the D matrix becomes (Timoshenko and Goodier, 1951), for unit thicknes EQ-») P= Tend - 2) 1 ° (264) 2.12. PLANE ELEMENT INERTIA MATRIX When inertia is significant (2.54) are supplemented by forces — p(@*u/ét*) and ~ p(g? jdt?) respectively, where pis the mass ofthe element per unit volume, For an clement of unit thickness this leads, in exactly the same way asin (2.14), to the element mass matrix which has terms given by Mot,=0 | [Nai deay 265) and hence to an eigenvalue equation the same as (2.19)30 213. AXISYMMETRIC STRESS AND STRAIN Solids of revolution subjected to axisymmetric loading possess only two independent components of displacement and can be analysed as if they were two-dimensional. For example, Figure 26(a) shows a thick tube subjected to radial pressure p and axial pressure q. Only atypical radial cross-section need be analysed and is subdivided into rectangular elements inthe figure. The cylindrical coordinate system, Figure 2.6b), is the most convenient and when itis used the clement stiffness equation equivalent to (2.58) is wor f f forpnearance 266) which, when integrated over one radian, becomes kat {[orpararae where the strain-displacement relations are now (Timoshenko and Goodier, 1951) oa ° LD, me 1 emde or a Figure 26 Cylinder under axial and radial pressure31 and as usual B= AN, where for elements of rectangular cross-section N could again be defined by (2.57) and (2.62). The stress-strain matrix must be redefined as E-» Tt (2.68) 2.14 THREE-DIMENSIONAL STRESS AND STRAIN When (2.48), (2.49) and (2.50) are extended to the three space-displacement variables u,v, , three simultaneous partial differential equations equivalent to (2.54) result. Discretisation proceeds as usual, and again the familiar element stiffness properties are derived as on [{ farppessyae (2.69) where the full strain-displacement relations are (Timoshenko and Goodier, 1951) a &) |x 9 0 @ 5, ° 5 a oo Zhu a = ° 2.70) 2 2 9 {lv tol lay oe aa ml | ap ay a a te) [5p 9 Fy or exAe and as before B= AN. For example, the rectangular brick-shaped element shown in Figure 2.7, which has eight corner nodes, would have shape functions of the32 type (+ ‘and so on. The full N matrix would be N, 0 0 N=/0 N, 0 (2.72) 0 ON =ON, Na Ns Ng Ns Ne Ny, Ne) (2.73) leading to the assembly of B. The stress-strain matrix is in this case any where x PTH =I a3 215 PLATE BENDING ELEMENT The bending of a thin plate is governed by the differential equation DV‘w=q (2.75) where D is the flexural rigidity of the plate, w its deflection in the transverse (2) direction and q the applied transverse load. The flexural rigidity is given by Eh? P= (2.76) where h is the plate thickness. Solution of equation (2.75) directly, for example by Galerkin minimisation, appears to imply that for a fixed D, a thin plate’s deflection is unaffected by the value of Poisson's ratio, This isin fact true only for certain boundary conditions and in general the integration by parts in the Galerkin process will supply extra terms which are dependent on ¥. This is a case in which the energy approach provides a simpler formulation. The train energy stored in a piece of bent plateis (Timoshenko and Woinowsky- Krieger, 1959) v0 f3) or una l (f(s) +S) ees #10022)" hacay @78) Consider, for example, the rectangular element shown in Figure 2.8, Ifthere are assumed to be four degrees of freedom per node, namely ow ow ew Fh ad Bae then the appropriate element shape functions can be shown to be products of the A Bay Figure 28 Rectangular plate bending element4 ‘beam shape functions already described. That is, as usual w=Nw 79) where, if the freedoms are numbered (W,0,. 45, 9xy)ace would be 234s the first term in N Bax? +23) 5 (6° -3by? +2y°) (2.80) rd Defining (@>—3ax? +23), Qy =e —3by? +2y°) (@x—2axt+°), Or=palby—2by*+ 9) 1 1 (281) Pyess(Bax?— 20, Qs mfp (Bby? 29") 1 Ppa al? ax), the full list of shape functions becomes N=P,Q,, N= P20s, Ny=Pi0s, N= P20, woP.0,, rr) Ne=P20s, N= PQs, Ne= P20. Using the same energy formulation terminology as before, define a M, ry ofS 2 M}ep}ot o |] & bw (283) a or M=DAw It can readily be verified that equation (2.78) for strain energy can be written u =f [iawrDiamae dy (284)35 and that the stiffness matrix becomes M= [fppsea, which is again the familiar equivalent of equation (2.47), with B= AN. ‘A typical value in the element stifiness matrix is given by [fOnseny , 22N22N, | dN, a KMy= if NEN, | ONOIN, | 22ND, Ge be FGF BF TYG Gy? NON ON, ON, ae yt 1 Naeby Bay ew) and these integrals are performed using Gaussian quadrature in two dimensions in Chapter 4, For some boundary conditions, the terms involving Poisson's r will cancel out ‘As usual, typical values in the mass matrix are given by Mat = ph | [vn,exay 87 2.16 SUMMARY OF ELEMENT EQUATIONS FOR SOLIDS ‘The preceding sections have demonstrated the essential similarity of al problems in solid mechanics when formulated in terms of finite elements. The statement of ‘element properties isto be found in two expressions, namely the element stifness matrix kM= f B'DBad (clement) (2.88) and the element mass matrix MM= ofsrwa (clement) (2.89) These expressions then appear in the three main classes of problem which concern us in engineering practice, namely () Static equilibrium problems KMr =F (2.90) i) Eigenvalue problems (KM —o?MM)r=0 291) ii) Propagation problems KMr +. ume =Fq) (292) Equations (2.90) are simultaneous equations which can be solved for known forces F to give equilibrium displacements r. Equations (2.91) may be solved by various techniques (iteration, QR algorithm, etc; see Bathe and Wilson (1976) or Jennings (1977) te yield mode shapes r and natural frequencies « of elastic36 systems, while equations (2.92) can be solved by advancing step by step in time from a known initial condition or by transformation to the frequency domain. Later chapters in the book describe programs that enable the user to solve practical engineering problems which are governed by these three basic equations. Additional features, such as treatment of non-linearity, damping, etc, will be dealt with in these chapters as they arise. 2.17 FLOW OF FLUIDS: NAVIER-STOKES EQUATIONS We shall be concerned only with the equations governing the motion of viscous, incompressible fluids. These equations are widely developed elsewhere, e Schlichting (1960), Preserving an analogy with previous sections on two- dimensional solids, u and v now become velocities in the x and y directions respectively and pis the mass density as before, Also as before, F, and F, are body forces in the appropriate directions. Conservation of mass leads to PaO On 042 i+ Z00)=0 293) ‘but due to incompressibility this may be reduced to = as) Conservation of momentum leads to au du ou 20, , Oey (2a rude ota + (Mee) 6 ax * 8 ot Ox Y, y (295) 20 4 12 4 9) gy (200 4 tn o( tend vo® ans (Sea be where 04,0, and t,, are stress components as previously defined for solids. Introducing the simplest constitutive parameters x (the molecular viscosity), 2 (taken to be — 3,1) and p (the fluid pressure) the following form of the stress equations is reached: , (du, av éu ao, av 1 ap_ gu (oo ov) _ usta pay a (at tat) ae Proceeding as before, and for the moment assuming the same shape functions are applied to all variables, u = Nu, v = Ny, p = Np we have for a single clement, treating the terms u and v as constants 7 = Nuy and = Nvo for the purposes of integration e298 299) aN, oN, LON, HON BON TUF Ut Png atm gar tee 100) aN, a gQN, p LON, HON _dIN aN Lan, "Ray Tp ay? poe" pay ut Multiplying by the weighting functions and integrating as usual yields an an Le yan nae, PN yacdy+4 { (NX pax ff a axdy+ | [wos d: y+ [In pdxdy #( (nN uaeay—4 [ (NEN, £[yPBaacay£| n2Sadsay-0 and Q101)38 ON {[waON Lf flan [JraSvava + | fro racay +i [raree -E [Jn Boast [nace Integrating products by parts where necessary and = resulting contour integrals gives éN ON 1 6N [JostBansyes [fotasyast il [fem H([ON ON 4, ‘ON ON +p lle sxteoret [lay ay teoren ww i tm [nstecons [fro narod [bac aN ON aN ON +h [FR Racaree BR ardya ‘The set of equations is completed by the continuity condition [Jn(eNer2s)esayea ay Collecting terms in u, p and v respectively leads to an equilibrium equation (Taylor and Hughes, 1981) Cy Ce (2.1048) where y= Cu [xavay 2.1040) ox a= aN Cua [pew39 Cy =0 1 aN Cy: = | |=N<—dxdy * Sf yaad C=C Referring to Table 2.1 we now have many terms ofthe type (ON, /@x) which imply unsymmetrical structures for Cy. Thus special solution algorithms will be necessary. Computational details are let until Chapter 9. 218 SIMPLIFIED FLOW EQUATIONS In many practical instances it may not be necessary to solve the complete coupled. system described in the previous section. The pressure p can be eliminated from (2.98) and if vorticity a is defined as ou_ a0 a (2105) this results in a single equation: (2.106) 2.107) ‘an alternative coupled system involving y and w can be devised, given here for steady state condition ay ev, ae a (2.108) e(coete\eweet a a p\ ax * Oy?) “dy ox ox dy This clearly has the advantage that only two unknowns are involved rather than the previous three. However, the solution of (2.108) isstill a relatively complicated process and flow problems are sometimes solved via equation (2.106) alone, assuming that wand v can be approximated by some independent means. In this form, equation (2.106) is an example of the ‘diffusion-convection’ equation, the second order space derivatives corresponding toa ‘diffusion’ process and the first order ones to a ‘convection’ process. The equation arises in various areas of40 engincering, for example sediment transport and pollutant disposal (Smith, 1976, 1979). If there is no convection, the resulting equation is of the type do _ 4 (Bo, Po fee (far eae 1 2 at (eS 2109) which s the ‘heat conduction’ or‘iffusion’ equation well known in many areas of engineering. ‘tinal simplification isa reduction to steady state conditions, in which case Bo Fo ional ati ae tay? en leaving the familiar ‘Laplace’ equation. In the following sections, finite element formulations of these simplified flow equations are described, in order of increasing complexity. 2.19 SIMPLIFIED FLUID FLOW: STEADY STATE, ‘The form of Laplace's equation (2.110) which arisesin geomechanics, for example concerning groundwater flow in an aquifer (Muskat, 1937), is 246 ab bath ss Quit Rather =O au where @ is the uid potential and k,,k, are permeabilities in the x and y directions. The finite clement discretisation process reduces the differential ‘equation to a set of equilibrium type simultaneous equations of the form KP¢=Q (2.112) where KP is symmetrical and Q is a vector of net nodal inflows/outflows. Reference to Table 2.1 shows that typical terms in the matrix KP are of the form ANON) 5 ON ON) [JOS enBBay ain With the usual finite element discretisation g=Nob (2114) a convenient way of expressing the matrix KP in (2.112) is xe [[rxraxay us) where the property matrix K is analogous to the stress strain matrix D in solid mechanics; thus k, 0 «=| a | (2.116)at {assuming that the axes of the permeability tensor coincide with x and y). The T matrix is similar to the B matrix of solid mechanics and is given by amy aN. éx Ox ON, ON, ONs ONe dy “ay “Oy oy The similarity between (2.115) for a fluid and (2.63) for a solid enables the corresponding programs to look similar in spite of the governing differential ‘equations being quite different. This unity of treatment is utilised in describing the programming techniques in Chapter 3 Finally, it is worth noting that (2.115) can also be arrived at from energy considerations. The equivalent energy statement is that the integral Sse) Gi) Jaw ey shall be a minimum for all possible (x, )). Example solutions to steady state problems described by (2.111) are given in Chapter 7. Quin 2.20 SIMPLIFIED FLUID FLOW: TRANSIENT STATE Transient conditions must be analysed in many physical situations, for example in the case of Terzaghi ‘consolidation’ in soil mechanics or heat conduction. The ‘governing consolidation diffusion equation for excess pore pressure 1, takes the form ly Pity _ Bity Saat +O OF oe (2.119) where ¢,,¢, ate the coefficients of consolidation. Discretisation of the left hand side of (2.119) clearly follows that of (2.111) while the time derivative will be associated with a matrix of the ‘mass matrix’ type without the multiple p. Hence the discretised system is (2.120) This set of first order, ordinary differential equations can be solved by many ‘methods, the simplest of which discretise the time derivative by finite differences. The algorithms are described in Chapter 3 with example solutions in Chapter 8 221 SIMPLIFIED FLUID FLOW WITH ADVECTION If pollutants, sediments, tracers, etc., are transported by a laminar flow system they are at the same time translated or ‘advected’ by the flow and diffused withi4a it. The governing differential equation for the two-dimensional case is (Smith etal, 1973) 74, 20 yb ob a8 2.121) Cxgya tr ayF “Gx Gy” Ot where wand v are thefluid velocity components the and y directions (compare equation 2.106) The extra advection terms — u(2¢/2x) and — (26/03) compared with (2.119) lead, es shown in Table 21, to unsymmetric components ofthe stiffness’ matrix of the type r {{( atten) an x ay which must be added to the symmetric, diffusion components given in (2113) ‘When this has been done, equilibrium equations like (2.112) or transient equa- tions like (2.120) are regained. ‘Mathematically, equation (2.121) is a differential equation which is not self- adjoint (Bere, 1962), due to the presence of the first order spatial derivatives From a finite element point of view, equations which are not self-adjoint will always lead to unsymmetrical stiffness matrices. "A second consequence of non-self-adjoint equations is that there is no energy formulation equivalent to (2.118) It is clearly a benefit of the Galerkin approach that it can be used for all types of equation and is not restricted to self-adjoint systems. ‘Equation (2.121) can be rendered self-adjoint by using the transformation (2.123) but thisis nor recommended unless wand v are small compared with c, and ¢, a shown by Smith etal, (1973) Equation (2.121) and the use of (2.123) are described in Chapter 8. 222 FURTHER COUPLED EQUATIONS: BIOT CONSOLIDATION ‘Thus far in this chapter, analyses of solids and fluids have been considered separately. However, Biot formulated the theory of coupled solid-fluid interac. tion which finds applications in soil mechanics (Smith and Hobbs, 1976). Thesoil skeleton is treated as a porous elastic solid and the laminar porefluid is coupled to the solid by the conditions of compressibility and of continuity. Thus Biot’s governing equation is given by KL, Oe yp Ole 4 g Oe selhas thar th ae? | where K’ is the soil bulk modulus and p is the mean total stress. (2.124)B For two-dimensional equilibrium in the absence of body forces, the gradient of fective stress from (2.48) must be augmented by the gradients of the fuid reSSUre thy a follows: Stay 5 Bile 0 y * Ox ens Fey 205 , Oty éx © dy * dy The constitutive laws are those previously defined for the solid and fluid respectively; hence in plane strain «, 1 o | ie Eun) | ¥ 2 acer Mi= 2) | venue nce (een 1-2v te 0 0 Al he Gl Ake 0) faua/ex 1 heads ele om where q, amd q, are the volumetric flow rates per unit area into and out of the element and 7, is the unit weight of water. The solid strain-displacement relations are still given by (2.50), and the final condition is that for full saturation and, in this case incompressibility, outflow from an element of soil equals the reduction in volume of the element. Hence d (du, av $ (H+) 2.128) and from (2.127) the third differential equation is given by ik, O71 @u, d (du av qe Ot te Op tae (+3)=° 1) Asusual ina displacement method, aand are eliminated in terms ofu and vso that the final coupled variables are u,v and u,. These are now discretised in the normal way: (2130) In practice, it may be preferable to use a higher order of discretisation for wand4 » compared with u, but, for the present, the same shape functions are used to describe all three variables. ‘When discretisation and the Galerkin process are completed, (2.125) and (2.128) lead to the pair of equilibrium and continuity equations: KMr + Cu, =F at rr where, for a four-noded element, cot KPu, =0 4 has ref} and wed 2.132) KM and KP are the elastic and fluid ‘stfiness’ matrices and C is a rectangular coupling matrix consisting of terms of the form [Jrctteox ay 2133) Fis the external loading vector. Equations (2.131) must be integrated in time by some method such as finite dilerences and this is described further in Chapter 3. Examples of such solutions in practice are given in Chapter 9. 2.23 CONCLUSIONS When viewed from a finite clement standpoint, all equilibrium problems, whether involving solids or fluids, take the same form, namely KMr=F or (2134) KP6=Q Forsimple uncoupled problems the solid KM and fluid KP matrices have similar symmetrical structures and computer programs to construct them will be similar. However, for coupled problems such as are described by the Navier-Stokes equations, KM (or KP) is unsymmetrical and appropriate alternative software will be necessary. Tn the same way, eigenvalue, propagation and transient problems all involve the mass matrix MM (or a simple multiple of it, PM), Therefore, coding of these45 different types of solution can be expected to contain sections common toall three problems. So far, single elements have been considered in the discretisation process, and only the simplest line and rectangular elements have been described. The next chapter is mainly devoted to a description of programming strategy, but before this, the finice element concept is extended to embrace meshes of interlinked elements and elements of general shape. 2.24 REFERENCES Bathe, K. J, and Wilson, E.L. (1976) Numerical Methods in Finite Element Analysis. Prentice-Hall, Englewood Clifs, New Jersey. Berg, P. N,(1962) Calculus of variations. In Handbook of Engineering Mechanics, Chapter 16, ed. W. Flugge, McGraw-Hill, New York Connor, J.J, and Brebbia, C.A."(1976) Finite Element Techniques for Fluid Flow. ‘Newnes-Butterworth, London Cook, RD. (1981) Concepts and Applications of Finite Element Analysis, 2nd edition, Wiley, New York, Crandall, 8, H. (1956) Engineering Analysis. McGraw-Hill, New York. Finlayson, B.A. (1972) The Method of Weighted Residuals and Variational Principles. ‘Academic Press, New York. Horne, M. R.,and Merchant, W. (1965) The Stability of Frames. Pergamon Press, Oxford, Jennings, A. (1977) Matrix Computation for Engineers and Scientists. John Wiley, London, Leckie, FA, and Lindberg, G. M. (1963) The effect of lumped parameters on beam frequencies, The Aeronautical Quarterly, 14, 234 Livesley, R. K. (1975) Matrix methods of structural analysis, Pergamon Press, Lundquist, E.E., and Kroll, W. D. (1944) NACA Report ARR, No. 4824 Muskat, M. (1937) The Flow of Homogeneous Fluids through Porous Media. McGraw-Hill, New York. Rao, S. S (1982) The Finite Element Method in Engineering, Pergamon Press, Oxford, Schlichting, H. (1960) Boundary Layer Theory. McGraw-Hill, New York Smith, I. M, (1976) Integration in time of diffusion and diffusion-convection equations. In Finite Elements in Water Resources ed. W.G. Gray, G. F. Pinder and C.A. Brebbia, pp. 1.3-1.20. Pentech Press ‘Smith, M. (1979) The difusion-convection equation. In A Survey of Numerical Methods for Partial Differential Equations, ed. 1. Gladwell and R. Wait, pp. 195-211. Oxford University Press Smith, L.M., Farraday, R.V., and O'Connor, B. A. (1973) Rayleigh-Ritz and Galerkin finite elements for diflusion-convection problems. Water Resources Research, 9, No.3, 593. ‘Smith, I. M,, and Hobbs, R. (1976) Biot analysis of consolidation beneath embankments, Géotechnique, 26, No. 1, 149. Strang, G., and Fix, G, J.(1973) An Analysis of he Finite Element Method. Prentice-Hall, Englewood Clifls, New Jersey Szabo, B. A. and Lee, G. C. (1969) Derivation of stiffness matrices for problems in plane clasticity by the Galerkin method. Im. J. Num. Meth. Eng. 1, 301 Taig, I. C. (1961) Structural analysis by the matrix displacement method, English Electric ‘Aviation Report No.SO17. ‘Taylor, C, and Hughes, T. G. (1981) Finite Element Programming of the Navier-Stokes Equation. Pineridge Press Ltd, Swansea,46 Timoshenko, & P, and Goodier, J.N. (1951) Theory of Elasticity. McGraw-Hill, New York. ‘Timoshenko, S.P,, and Woinowsky-Krieger, 8. (1959) Theory of Plates and Shells. ‘McGraw-Hill, New York Zienkiewiez, O.C(1977) The Finite Element Method in Engineering Science, 3rd edition. ‘McGraw-Hill, London.CHAPTER 3 Programming Finite Element Computations 3.0 INTRODUCTION In Chapter 2, the finite element discretisation process was described, whereby partial differential equations can be replaced by matrix equations which take the form of linear and nonlinear algebraic equations, cigenvalue equations or ordinary differential equations inthe time variable. The present chapter describes hhow programs can be constructed in order to formulate and solve these kinds of equations. Before this, two additional features must be introduced, First, we have so far dealt only with the simplest shapes of elements, namely lines and rectangles. Obviously if differential equations are to be solved over regions of general shape, elements must be allowed to assume general shapes as well, This is accomplished by introducing general triangular and quadrilateral elements together with the concept of a coordinate system local to the element. Second, we have so far considered only a single element, whereas useful solutions will normally be obtained by many elements, usually hundreds or thousands in practice, joined together at the nodes. In addition, various types of boundary conditions may be prescribed which constrain the solution in some way. Local coordinate systems, element assembly and incorporation of boundary Conditions are all explained in the sections that follow. 34 LOCAL COORDINATES FOR QUADRILATERAL ELEMENTS Figure 3.1 shows two types of four-noded quadrilateral clements. The shape functions for the rectangle were shown to be given by equations (257), namely Nj, =(1—x/a)(1 ~ y/b) and so on. If it is attempted to construct similar shape functions in the ‘global’ coordinates (x,y) for the general quadrilateral, very complex algebraic expressions will result, which are at the very least tedious to check (Irons and Ahmad, 1980) Instead it is better to work in a local coordinate system as shown in Figure 3.2, originally proposed by Taig (1961). The general point P(é,n) within the a7y y (0,6) $5 : a 4 Toor x x Figure 3.1. (a) Plane rectangular element. (b) Plane general quadrilateral clement Figure 32 Local coordinate system for quadrilateral elements ‘quadrilateral is located at the intersection of two lines which cut opposite sides of the quadrilateral in equal proportions. For reasons associated with subsequent numerical integrations it proves to be convenient to ‘normalise" the coordinates so that side 1 2has €= 1, side 3 4 has = +I, side 1 4 has = — 1 and side 23 has n= + 1. In this system the intersection of the bisectors of opposite sides of the quadrilateral is the point (0,0), while the corners 1, 2, 3,4 are (1, ~ 1) (1.0, (1,0, (1,1) respectively. ‘When this choice is adopted, the shape functions for a four-noded quadrilateral with corner nodes take the simple form: Ny=KI-O(1~n) 41~O+n) (+ +n) =h1+O(l—n) and these can be used to describe the variation of unknowns such as displacement or fluid potential in an element as before. GB.)4” Under special circumstances the same shape functions can also be used to specify the relation between the global (x,y) and local (f,n) coordinate systems. If this is so the element is of a type called ‘isoparametric’ (Ergatoudis et a, 1968; Zienkiewice et al, 1969); the four-node quadrilateral is an example. The coordinate transformation is therefore XEN yx, + Nata + NgXy + Naty =Nx Y= Nay + Nav + Nays + Nave =Ny (3.2) where the N are given by (3.1) In the previous chapter it was shown that element properties involve not only N but also their derivatives with respect to the global coordinates (x, )) which appear in matrices such as B and T. Further, products of these quantities need to be integrated over the element area or volume. Derivatives are easily converted from one coordinate system to the other by means of the chain rule of partial differentiation, best expressed in matrix form by (a) [oe (a ag] ae jax, - J : jaf ~ ax 2 a (én) én ay, e fe) fa ax. og a. le (3.4) Gy an, where J is the Jacobian matrix. The determinant of this matrix, det J}, must also be evaluated because it is used in the transformed integrals as follows: fora f° fi det || dg dy G5) J BS Figure 33 (a) Degenerate quadrilateral. (b) Unacceptable ‘quadrilateralTable 3.1 Coordinates and weights in Gaus- sian quadrature formulae . ty wy 2 + 7 1 sh oo Under certain circumstances, for example that shown in Figure 3.3(b), the Jacobian becomes indeterminate. When using quadrilateral elements, reflex interior angles should be avoided. 32 NUMERICAL INTEGRATION FOR QUADRILATERALS Although some integrals of this type could be evaluated analytically, it is impractical for complicated functions, particularly in the general case when (&,7) become curvilinear. In practice (3:5) are evaluated numerically, using Gaussian quadrature over quadrilateral regions (Irons, 19662, 1966b). The quadrature rules in two dimensions are all of the form j fi senvecar= x Smo fon) G6) where w, and w, are weighting coefficients and ,, n; ate coordinate positions within the element. These values for n equal to 1,2 and 3 are shown in Table 3.1 and complete tables are available in other sources, e.g. Kopal (1961). The table assumes that the range of integration is + 1; hence the reason for normalising the Tocal coordinate system in this way. ‘The approximate equality in (36) is exact for cubie functions when n= 2 and for quintics when n = 3. Usually one attempts to perform integrations over finite clements exactly, but in special circumstances (Zienkiewicz et al, 1971) ‘redueed’ integration whereby integrals are evaluated approximately can improve the quality of solutions. 33 LOCAL COORDINATES FOR TRIANGULAR ELEMENTS Figure 3.4 shows how a general triangular clement can be mapped into a right- ‘angled isosceles triangle. This approach is identical to using area coordinates (Zienkiewiez ex al, 1971) in which any point within the triangle can be referenced using three local coordinates (Ly,3, 3). Clearly for a plane region, only (wost (107 by Figure 34 (a) General triangular element. (b) Element mapped. in local coordinates independent coordinates are necessary; hence the third coordinate isa function of the other two: Ly=1-Lj-ky @7 However, it often leads to more elegant formulations algebraically if all three coordinates are retained. For example, the shape functions for a three-noded triangular element (constant strai 3.8) and the isoparametric property is that Nx, + Nake + Nay 3s Y= Nays + Nays + Nays a Equations (3.3) and (3.4 from the previous paragraph still apply regarding the Jacobian matrix but equation (3.5) must be modified for triangles to give Jao ff 34 NUMERICAL INTEGRATION FOR TRIANGLES “det|S|dLyaL, G10) ‘Numerical integration over triangular regions s similar to that for quadrilaterals except that the sampling points do not occur in such ‘convenient’ positions. ‘The quadrature rules take the general form 1 pits A f f I (lysLa)dba dL 4 WSL Li) G.I) where w, is the weighting coefficient corresponding to the sampling point (Li, L5) and n represents the number of sampling points. Typical values of the weights and sampling points are given in Table 3.2.2 Table 32 Coordinates and weights for integr- ation over triangular areas ‘As with quadrilaterals, numerical integration can be exact for certain polynomials. For example, in Table 3.2, the one-point ruleis exact for integration of first order polynomials and the three-point rule is exact for polynomials of second order. Reduced integration can again be beneficial in some instances. ‘Computer formulations involving local coordinates, transformation of co ordinates and numerical integration are described in subsequent paragraphs, 35 ELEMENT ASSEMBLY Properties of elements in isolation have been shown to be giv equations, for example the equilibrium equation KP¢=Q G.12) describing steady laminar fluid flow. In Figure 3.5 is shown a small mesh containing three quadrilateral elements, all of which have properties defined by G.12). The next problem isto assemble the elements and so derive the properties of the three-element system, Each element possesses node numbers, shown in parentheses, which follow the scheme in Figures 3.1 to 3.3, namely numbering clockwise from the lower let hand corner. Since there is only one unknown at by matrix 7 . Figure 35. Mesh of quadilater! elements33 every node, the fluid potential, each individual element equation can be written KPia KPia KPi3 KPia) (oi) (0: KP2, KPz2 KP23 KP2al Joa|_ |Q2 KP3, KP32 KPs3 KPy.J )@s(~ Qs a KPa, KPs2 KPsx KPsa} {os} [Qs However, in the mesh numbering system, not in parentheses, mesh node 4 corresponds to element node (I) in element I and to element node 2) in element 3. The total number of equations for the mesh is 8 and, within this system, term KP, from element | and term KP, from element 3 would be added together and Would appear in location 4,4 and so on. The total system matrix for Figure 3.5 is given in Table 33, where the superscripts refer to element numbers. ‘This total system matzix is symmetrical provided its constituent matrices are symmetrical. The matrix also possesses the useful property of “bandedness’ Which means that the terms are concentrated around the ‘Teading diagonal which stretches from the upper left to the lower right ofthe table. In fact no term in any row can be more than four locations removed from the leading diagonal so the system is said to have a ‘semi-bandwidth’ IW of 4. This can be obtained by inspection from Figure 3.5 by subtracting the lowest from the highest freedom. number in each element. Complicated meshes have variable bandwidths and useful computer programs make use of banding when storing the system matrices The importance of efficient mesh numbering is illustrated for a mesh of line elements in Figure 3.6 where the scheme in parentheses has [W = 13 compared with the other scheme with 11V- Ifsystem symmetry exists it should also be taken into account. For example, if the system in Table 3.3 is symmetrical, there are only 30 unique components (the leading diagonal terms plus @ maximum of four terms to the left or right of the diagonal in each row). Often, with a slight decrease of efficiency, the symmetrical half of a band matrix is stored as a rectangular atray with a size equal to the ‘Table 3.3. System stiffness matrix for mesh in Figure 3.5, Superscripts indicate element ‘numbers KPL,_ KPhy 0 KPA, KPLa o 0 0 KP, KP\y+KP}; KP3; KP}, KPIs +KP3, KP}, 0 0 ° KP3.KP33 ° KPi, KPa 00 KPh, PL; 0. KPI, +KP2y KPLa+KPly 0. KPI, KPLy KPi, KP} +KPt, KP}, KPi,+KP32 KPi4+KPi, KPi4 KP3. KPia +KPRy 0 KPi, KPis oO KRY KPi, 0 oO oO 0 oO KP}: KPhy 0 KP, KPis o 0 0 KP, KPLS 0 KPL, KP,Figure 36 Alternative mesh ‘numbering schemes ‘number of system equations times the semi-bandwidth plus 1. In this case zeros are filled into the extra locations in the first few or last few rows depending. fon whether the lower or upper “half” of the matrix is stored. Special storage schemes, for example ‘skyline’ techniques (Bathe and Wilson, 1976), are also considered later in the chapter (see Figure 3.17) and are important in cases where bandwidths vary greatly Later in this chapter, procedures are described whereby system matrices like that in Table 3.3 can be automatically assembled in band form, with or without the ‘skyline’ option, from the constituent element matrices, 36 INCORPORATION OF BOUNDARY CONDITIONS, Eigenvalues of frecly floating clements or meshes are sometimes required but normally in cigenvalue problems and always in equilibrium and propagation problems additional boundary information has to be supplied before solutions ‘can be obtained. For example, the system matrix defined in Table 3.3 is singulat and this set of equations has no solution The simplest type of boundary condition occurs when the dependent variable in the solution is known to be Zero at various points in the region (and hence nodes in the finite element mesh). When this occurs, the equation components associated with these nodes are not required in the solution and information is given (o the assembly routine which prevents these components from ever being assembled into the final system. Thus only the non-zero nodal values are solved for. ‘A variation of this condition occurs when the dependent variable has known, ‘but non-zero, values at various locations. Although an elimination procedure could be devised, the way this condition is handled in practice is by adding a35 ‘targe’ number, say 10?.to the leading diagonal ofthe ‘stiffness’ matrix in the row in which the prescribed value is required. The term in the same row of the Fight hand side vector is then set to the prescribed value multiplied by the augmented “atffness’ coefficient, For example, suppose the value ofthe fluid head at node Sin Figure 3.5 is known to be S70units. The unconstrained set of equations (Table 3.3) would be assembled and term (5,5) augmented by adding 10", In the subsequent solution there would be an equation (Kg.s-+ 102} + small terms = $7.0 x (Ks,5 + 10") G.14) which would have the effect of making ¢s equal to 57.0. Clearly this procedure is only successful if indeed ‘small terms’ are stall relative to 10", ‘Boundary conditions can also involve gradients of the unknown in the forms og in 7° 15) ob moe .16) ob a 7, where 1 is the normal to the boundary and C,, C, are constants. To be specific, consider a solution of the diffusion-advection equation (2121) subject to boundary conditions (3.15), (3.16) and (3.17) respectively. When the second order terms c,(276/dx2) and c,(02p/0y) are integrated by parts, boundary integrals of the type Jenene G18) aise, where s isa length of boundary and I, the direction cosine of the normal. Clearly the case d¢/0n = 0 presents no difficulty since the contour integral (3.18) vanishes. However, (3.16) gives rise to an extra integral, which for the boundary element shown in Figure 3.7 is 5 f NTC, as 6.19) ) When ¢ is expanded as N@ we get an additional matrix oo00 (x,—x)}0 21 0 ; 3%) 6 o120 620) 000 0 Which must be added to the left hand side of the element equations.56 Figure 37 Boundary conditions involy- ing non-zero gradients of the unknown For boundary condition (3.17) the additional term is, ferrets 62) which is just a vector 0 Cxeadu=y) J0 —aT NN (3.22) 1 which would be added to the right hand side of the element equations. For a further discussion of boundary conditions see Smith (1979. In summary, boundary conditions of the type @ = 0 oF 04/0n=0 are the most common and are easily handled in finite element analyses, The cases = constant or @{dn = constant x ¢ are somewhat more complicated but can be appro- Priately treated. Examples of the use of these types of boundary specification are included in the examples chapters. 3.7 PROGRAMMING USING BUILDING BLOCKS The programs in subsequent chapters are constituted from 100 or so building blocks in the form of FORTRAN subroutines which perform the tasks of ‘computing and integrating the element matrices, assembling these into system ‘matrices and carrying out the appropriate equilibrium, eigenvalue or propag- ation calculations. 11 is anticipated that most users will elect to pre-compile all of the building blocks and to hold these permanently in a library on backing store (usually disc}. The library should then be automatically accessible to the calling programs by ‘means of a simple attach in the job description. A summary of the building blocks is listed in Appendix 4 where their actions and input/output parameters are described, A separation has been made into “black box’ routines (concerned with some matrix operations). whose mode of37 action the reader need not necessarily know in detail, and special purpose routines which ate the basis of finite element computations. These special purpose routines are listed in Appendix 5. The black box routines should be thought ofas an addition to the permanent library functions such as SIN or ABS, ‘and could well be replaced from a mathematical subroutine library such as NAG. 38 BLACK BOX ROUTINES ‘The first group of subroutines is concerned with standard matrix handling, and. the action of these should be self-explanatory. The routines, usually assuming two-dimensional matrices, are MATCOP (copies) MATRAN (transposes).- MATSUB (subtracts) MATMUL (multiplies) NULL (nulls) TWOBY? (inverts 2 x 2) MVMULT (multiplies MATADD (adds) ‘TREEX3 (inverts 3 x 3) matrix by vector) MATINV (inverts N x N) MSMULT (multiplies NULLS (3-d array) matrix by scatar) Similar routines associated with vectors are VECCOP VVMULT (cross-product) VECADD NULVEC The second batch of subroutines is concerned with the solution of linear algebraic equations (required in equilibrium and propagation problems). The subroutines have been split into reduction and forward/backward re-substitution phases. BANRED CHOLIN — GAUSBA BACK] CHOBKI —SOLVBA BACK? CHOBK2 Method Gauss, Choleski Gauss Equation Symmetrical Symmetrical _Unsymmetrical coefficients COMRED — SPARIN COMBK1 —SPARBI COMBK2 PARR? Gauss Choleski Symmetrical Symmetrical complex skyline In the majority of programs, the forward and backward substitution is enshrined in the single subroutines BACSUB, CHOBAC, COMBAC and SPABAC respectively Several subroutines are associated with eigenvalue and cigenvector determin- ation, for example58 Tridiagonalisation — TRIDIA BANDRD. Finds eigenvalues. © EVECTS BISECT Method QR factorisation Jacobi Equation coefficients | Symmetrical ‘Symmetrical, banded, full Nx N. upper triangle stored stored ‘The first of each pair tridiagonalises the matrix and the second extracts all of the eigenvalues. The eigenvectors re found only by EVECTS. Itshould be noted that these routines, although very robust and accurate, can be inefficient both in storage requirements and in run-time and should not be used for solving very large problems, in which in any caseit is unlikely that the full range of eigenmodes would be required. The various veetor iteration methods (Bathe and Wilson, 1976) should be resorted to in those cases. ‘One of the most effective of these is the Lanczos method, and routines LANCZ1 LANCZ2 are used to calculate the Lanezos vectors and eigenvectors of a matrix. Although simple matrix-by-vector multiplications can be accomplished by MVMULT, advantage is usually taken of banding of the matrix coefficients ‘whenever possible. To allow for this, four special matrix-by-veetor multiplication routines are provided: LINMUL BANMUL BANTML _LINMLS Matrix coefficients Symmetrical Symmetrical Unsymmetrical Symmetrical Storage of matrix Vector Lower triangle Full band Skyline Further information on when these routines should be used is given in Section 3.10. During eigenvalue computation a couple of special matrix operations are catered for by LBKBAN _ Lest whose action is described in Appendix 4. Ina teaching text such as this, elaborate input and output procedures are avoided. It is expected that users may pre-process their input data using independent programs and plot the output graphically, At present’ these operations are very machine-ependent. Therefore simple free format input to “TAPE Sis adopted in all programs, The results of calculations are usually vectors of arrays and two subroutines enable the output of these in a standard E format to “TAPE 6. They are called PRINTV* PRINTA respectively.9 In order to describe the action of the remaining special purpose subroutines, which are listed in full in Appendix 5, t is necessary first to consider the properties of individual finite elements and then the representation of continua from assemblages of these elements Static linear problems (including eigenprablems) are considered first. Thereafter modifications to programs to incorporate time dependence are added. 39 SPECIAL PURPOSE ROUTI ‘The job of these routines is to compute the clement matrix coefficients, for example the ‘stiffness’, to integrate these over the element area or volume and finally to assemble the element submatrices into the global system matrix or matrices. The black box routines for equation solution, eigenvalue determination and so on then take over to produce the final results 39.1 Element matrix calculation In the remainder of this chapter the notation adopted is that used in the subroutine listings in Appendix 5. Wherever possible, mnemonics are used so that local coordinate ¢ becomes XI in the subroutines and so on. 39.1.1 Plane strain (stress) analysis of elastic solids using quadrilateral elements ‘As an example of element matrix calculation, consider the computation of the element stiffness matrix for plane elasticity given by (2.63): kM= ff B'DBdxdy 623) In program terminology this becomes xat= { fase De «BEE dxdy (3.24) and its formation is described by the inner loop of the structure chart in Figure 3.8, Itis assumed for the moment that the element nodal coordinates (x,y) have been calculated and stored in the array COORD. For example, we would have for a four-node quadrilateral CooRD =|*? 2? (3.25) Os 4. ‘The shape functions N are held in array FUN, in terms of local coordinates, asFindthe element geometry (nodal coordinates, COORD) Nullthestitfooss matrixspace, KM Forallthe Gaussian integrating points, NGP Find the Gauss point coordinates and weighting factorsin array SAMP. Formthe element shape functions, FUN and dorivates, DERin local coordinates XI, ETA. Global coordinate ranstormation of derivatives toDERWV. Formthe strain displacement matrix, BEE “Multiply by he stess strain matrix, OBE. togive DBEE, Transpose BEE to BT andmultply into BEE togive BTDB, ‘Add this contribution tothe element stifness KM. ‘Agsomble the current KMinto the global system matrix. Figure 38 Structure chart of element matrix assembly specified in (3.1) by jd—xna—eTA))* 41X11 + ETA) 4+X0(1+ ETA) 40+ Xb(1— ETA) ‘The B matrix contains derivatives of the shape functions and these are easily computed in the local coordinate system as 3.26) ora per=!{-(-ETA) -(1+ETA) (1+ETA) (1—ETA) -(=-xX) (=X) +X) 0 +X 4] ‘The information in (3.26) and (3.27) for a four-node quadrilateral is formed by the subroutine | 27, FORMLN for the specific Gaussian integration points (XI, ETA), held in the array SAMP where I and J run from I to NGP, the number of Gauss points specified in each direction. Figure 3.9(a) shows the typical layout und ordering for two-point Gaussian integration, In all-cases SAMP is formed by the subroutine Gauss where NGP can take any value from 1 t0 7. The derivatives DER must then be converted into their counterparts in the (x,y) coordinate system, DERTV, by means of the Jacobian matrix transform= ation (3.3) or (3.4). From the isoparametric property, {F}-coorprerunt (328) and since the Jacobian matrix is given by zig % a 3.29) eoy a on On itis clear that (3.29) can be obtained by dillerentiating (3.28) with respect to the local coordinates. In this way JAC = DER+COORD (330) Inorder to compute DERIV we must invert JAC to give JACI using TWOBY2 1 u it ' re bs lz, tae Ba++ | NGPe2 NIP=3 u (0) (b) Figure 39 Integration schemes for (a) quadrilaterals and (b) trianglesa2 in this two-dimensional case and finally carry out the multipli DERIV =JACI+DER Gat) Thus the sequence of operations CALL FORMLN (DER, IDER, FUN, SAMP, ISAMP, I, J) CALL MATMUL (DER, IDER, COORD, ICOORD, JAC, UAC, IT, NOD, IT) CALL TWOBY? JAC, DAC, JACI, JACI, DET) 632) CALL MATMUL (JACI, ACI, DER, IDER, DERIV, IDERIY, IT, IT, NOD) where NOD is the number of element nodes four in this case) and IT the number of coordinate dimensions (two in this case), will be found in programs for plane elasticity using four-noded quadrilateral elements. After these operations have been performed, the derivatives of the element shape functions with respect to (x,y) are held in DERIV while DET is the determinant of the Jacobian matrix, required later for the purposes of numerical integration. Following the conven tion adopted in Chapter 1, IDER is the working size of array DER and so on. “The matrix BEE in (3.24) can now be assembled asit consists of components of DERIV. This assembly is performed by the subroutine FORMB for plane problems. Thus the strain-displacement relations are EPS = BEE+ELD 633) where in the case of a foursnode quadrilateral ELD = {uy 05 ty Py My 05 ta 04) 34) “The variables wand vare simply the nodal displacements inthe x and y directions respectively assuming the nodal ordering of Figure 35 ‘The components of the integral of BEE*DEE* BEE, evaluated at the Gauss points given by all combinations of | and J,can now be computed by transposing BEE (o give BT, by forming the stress-strain matrix using the subroutine FMDEPS (for plane strain) where (335)6 and by carrying out the multiplications BTDB = BT+DEE*BEE 336) AA plane stress analysis would be obtained by simply replacing FMDEPS by FMDSIG. The integral is evaluated numerically by KM=DET+ ¥ WieWy«BTDB, 337) where W, and W, are the Gaussian weighting coeflicients held in array SAMP. ‘As soon as the element matrix is formed from (3.37) itis assembled into the global system matrix (or matrices) by special subroutines described later in this chapter. Since the assembly process is common to all problems, modifications to the clement matrix calculation for different situations will first be described. 39.1.2 Plane strain (stress) analysis of elastic solids using triangular elements ‘The previous section showed how the stiffness matrix of a typical four-node quadrilateral could be built up. In order to use triangular elements, very few alterations are required. For example, for a six-node triangular element, the values of the shape functions and their derivatives with respect to local coordinates at a particular location (Ly,L2,Ls) are formed by the subroutine FMTRI6 “This delivers the shape functions @L,~1)L, )* 4L,L; _} @b,—DL: Fue} Oa) 6.38) Ly — 1)Ly 4L4L, and their derivatives with respect to L, and Ly aFUNT peR=|_°" @FUNT ay [eho 4, 0 =4L; -@Ly-1) 4(ly-L) 0 4Ly Ly-1) 4(Ly-L,) -@Ly-1) AL, (39) The nodal numbering and the order in which the integration points are sampled for a typical three-point scheme are shown in Figure 3.9(b). ‘64 For integration over triangles, the sampling points in local coordinates (Lj,L,) are held in the array SAMP and the corresponding weighting cocficients in the vector WT. Both of these items are provided by the subroutine NUMINT ‘The version ofthis subroutine described in Appendix $ allows the total number of integrating points (NIP) to take the value I, 3, 4 6, 7, 12 or 16. The coding should be referred to in order to determine the sequence in which the integrating points are sampled for NIP > 3. The sequence of operations CALL FMTRI6 (DER, IDER, FUN, SAMP, ISAMP, I) CALL MATMUL (DER, IDER, COORD, ICOORD, JAC, JAC, IT, NOD, IT) G40) CALL TWOBY? (JAC, AC JACI, JACI, DET) CALL MATMUL (JACI, IACI, DER, IDER, DERIY, IDERIV, IT, IT, NOD) (here NOD now equals 6) places the required derivatives with respect to (x, y)in DERIV and finds the Jacobian determinant DET. Finally, numerical integration is performed by KM=05«DET+ }. W+BTDB, Gal) the factor of 0.5 being required because the weights add up to unity whereas the ‘tea of the triangle in local coordinates only equals one half ‘A higher order triangular element with fifteen nodes is also considered in Chapter 5. The shape function and derivatives for this element are provided by routine FMTRI5; otherwise the sequence of operations is virtually identical to those described above for the six-node element. 39.1.3 Axisymmetric strain of elastic solids The strain-displacement relations can again be written by (3.33) but in this case BEE must be formed by the subroutine FMBRAD where the cylindrical coordinates (r,2) replace their plane strain counterparts (4.9), The stress-strain matrix is still given by an expression similar to (3.35) but the 4x 4 DEE matrix is formed by FMDRAD In this case the integrated element stiffness is (2.66), namely KM -| JorsDerenEEsrara: 642)65 Considering a four-node element, the isoparametric property gives, r=FUN 1 "4 = FUN(K)}+COORDEK) =SUM 43) Hence we have KM=SUM«DET# YY) WieW;BTD3,, 44) By comparison with (3.37) it may be seen that when evaluated numerically the algorithms for axisymmetric and plane stiffness formation will be essentially the same, despite the fact that they are algebraically quite different. This is very significant from the points of view of programming effort and of program flexibility. However (3.44) now involves numerical evaluation of integrals involving l/r which do not have simple polynomial representations. Therefore, unlike planar problems, it will be impossible to evaluate (3.44) exactly by numerical means, especially as r (ie. SUM) approaches zero, Provided integration points do not lie onther = Oaxis, however, reasonable results are usually achieved using a similar order of quadrature to planar analysis. 39.14 Plane steady laminar fluid flow It was shown in (2.115) that a fluid element has a ‘stiffness’ defined by xo { frrxrasay 45) which becomes in program terminology x= [fosuvsexayepenivessy a and the similarity to (3.24) is obvious. The matrix DERIV simply contains the derivatives of the element shape functions with respect to (x,3) which were previously needed in the analysis of solids and are formed by the sequence (3.32) while KAY contains the permeability properties of the clement in the form PX 0 [* 2] oo In the computations, DTKD = DERIV'eKAY+DERIV G48) KAY6 is formed by appropriate matrix multiplications and the final matrix summation for a quadrilateral element is KP=DETe SS. WeWysDTKD,, 49) By comparison with (3.37) it will be seen that these physically very different problems are likely to require similar solution algorithms. 3.9.15. Mass matrix formation ‘The mass or inertia matrix was shown in Chapter 2, e.g, (2.65), to take the general form Mi=p | [NTNaxay 3.50) where N are just the shape functions. In the case of plane Muid flow, since there is only one degree of freedom per node, the ‘mass’ matrix is particularly simple in program terminology, namely MM ~ RHO» | [FUNT+FUNaxdy G51) By defining the product of the shape functions FTF = FUN'*FUN 52) we have MM =RHO*DET+ SY) W.eWyeFTF, G53) where RHO is the mass density. In the case of plane stress or strain of solids, because of the arrangement of the displacement vector in (3.34) it is convenient to use a special subroutine ECMAT to form the terms of the mass matrix as ECM before integration. Thereafter Mat = RHO. { [eM axay =RHO+DET+ 5.) Wy#W,+ECM,, G54) When ‘iumped’ mass approximations are used MM becomes a diagonal matrix. For a four-noded quadrilateral (NOD = 4), for example, MM =(RHO+AREA/NOD)«I 3.55) where AREA is the clement area and I the unit matrix. For higher order elements, however, all nodes may not receive equal weighting.3.9.1.6 Higher order quadrilateral elements 67 To emphasise the ease with which element types can be interchanged in programs, consider the next member of the isoparametric quadrilateral group, namely the ‘quadratic’ quadritateral with midside nodes shown in Figure 3.10 The same local coordinate system is retained and the coordinate matrix becomes % coorp = | ** xs X ‘The shape functions are now 4(1—XD(1 — ETA)(—XI- ETA —1) 41 =X1)(1 ETA?) 4= X11 + ETA)(—XI + ETA ~1) run =] {0-XP)0+ ETA) (1 +X(1 + ETA)(X1+ ETA~1) 4 +X1)(| — ETA?) (1 +X0(1 — BTA)(XI- ETA ~1) 4q_—X1-)(1— ETA) yi Ya vs Ye vs Ye yy vr 3.56) 57) which, together with their derivatives with respect to local coordinates, DER, are Figure 3.10 General quadratic quadsilateral ‘element68 formed by the subroutine FMQUAD The sequence of operations: CALL FMQUAD (DER, IDER, FUN,SAMP, ISAMP, I,J) CALL MATMUL (DER, IDER, COORD, ICOORD, JAG, JAC, IT, NOD, IT) CALL TWOBY?2 (JAC, JAC, JACI, JACI, DET) 58) CALL MATMUL (ACI, ACI, DER, IDER, DERIV, IDERIV, IT, IT, NOD) (where NOD now equals 8) places the required derivatives with respect to (x,y) in DERIV and finds the Jacobian determinant DET. ‘Another plane clement used in the programs later in this book is the Lagrangian nine-node element, This element uses ‘complete’ polynomial inter- polation in each direction, but requires a ninth node at its centre (Figure 3.11). ‘The shape Functions for this element are given by $(X1)(X1 = 1)(ETAV(ETA — 1) t = }XI(XT— (ETA + (ETA = 1) HOXI(XI— (ETAMETA +10) —HOXI + 1X1 — 1)(BTAV(ETA +1) FUN=} HXI\XI+ N(ETAWETA +1) G59) = HXI)(XI-+ ETA + (ETA 1) FOXX + I(ETA\ETA ~ 1) = HOX1 + 1X1 — 1(ETANETA = 1) (X14 (XI N(ETA + 1(ETA ~1) Figure 3.11 The Lagrangian nine-node ‘element9 and are formed together with their derivatives with respect to local coordinates by the subroutine FMLAG9 A comparison of (3.58), (3.40) and (3.32) indicates that programs using diferent element types will be almost identical, although operating on different sizes of arrays. Indeed, the reader could readily devise a completely general two- dimensional program in which element type selection was accomplished by data (ie, FORMLN, FMQUAD, FMLAG9, FMTRIG, ete, could be merged into @ single routine) 39.17 Three-dimensional cubic elements ‘As was the case with changes of plane element types, changes of element dimensions are readily made. For example, the eight-node brick element in Figure 3.12 is the three-dimensional extension of the four-noded quadrilateral Using the local coordinate system (7,0) the coordinate matrix is HOM Hi Sab Xs Ya fa coorp=|*+ % 7 (2.602) Bs Ys 25 Xe Ye 4 Xe Xe Ye Fe Figure 3.2 General linear brick element0 ‘The shape functions are {1 —X1)(1 = ETA) = ZETA) | 41 — X(t — ETA\(I + ZETA) [(1 + XI(1 = ETA) + ZETA) 4(.+XI01— ETAV(1 ~ZETAD | —X1)(1 + ETA)(I -ZETA) | 41 —XD(. + ETA)(I + ZETA) AU+XD(.+ ETA) + ZETA) 4 + XIQ+ ETA)(1 ~ZETA) (3.606) which together with their derivatives with respect to local coordinates are formed by the subroutine FMLIN3 The sequence of operations: CALL FMLIN3 (DER, IDER, FUN, SAMP, ISAMP.1.J,K) CALL MATMUL (DER, IDER, COORD, ICOORD, JAC, JAG, IT, NOD, IT) CALL TREEX3 (JAC, JAC, JACI, IACI, DET) (361) CALL MATMUL (JACI, JACI, DER, IDER, DERIV, IDERIV, IT, IT, NOD) (where NOD now equals 8 and IT now equals 3) results in DERIV, the required gradients with respect to (x,y,2) and the Jacobian determinant DET. rer order brick element with 20 nodes (Figure 3.13) is also used in x : Bie 4 3 18 2 19 1 2 7 Figure 313 A 20-node brick elementn programs later in the book. The shape functions for this element are provided by the subroutine FMQUA3 and are not quoted here, but see the coding in Appendix 5. For the three-dimensional elastic solid the element stiffness is given by ror= [[forsoresass dca a where BEE and DEE must now be formed by the subroutines FORMB3 \ FORMD3 respectively. The final summation becomes KM=DET+ )" SY. WytWye WesBTDByy x (3.63) where W,, W,, Wx are as usual Gaussian multipliers obtained from SAMP. For three-dimensional steady laminar fluid flow the element ‘stiffness’ is KP= il} |DERIV'-KAY«DERIV dxdy dz “ff ‘pal axes permeability tensor forkpaxayas cr) where KAY is the p PX 0 0 . Kay=|0 PY 0 3.65) 0 0 PZ Gaussian integration gives Nap NaF Nor KP=DET« SY, WieWeWaeDTKD is (3.66) Which is similar to (3.63). ‘The mass matrix for potential low is MM = RHO+ [front-FuNasayas =RHO+ [[[rreaxayas es Which is replaced by quadrature as MM = RHO+DET+ ¥° > WeWeWxtFTFy x (3.68)n In solid mechanics, a three-dimensional equivalent of ECMAT would be Fequited. This development is left to the reader. 39.18 Three-dimensional tetrahedron elements ‘An alternative clement for three-dimensional analysis is the tetrahedron, the simplest of which has four corner nodes. The local coordinate system makes use of volume coordinates as shown in Figure 3.14, For example, point P can be identified by four volume coordinates (Ly,L3,Ls, La) where 432 3.69) As is to be expected, one of these coordinates is redundant due to the identity Lytlythyt+ly=1 6.70) The shape functions for the constant strain tetrahedron are lk Ly FUN= 1)? an) L, z yf x 2 Volume of 1 tetrahedron = Vy Figure 3.14 A four-node tetrahedron element