PROGRAMMING THE FINITE ELEMENT METHOD SECOND EDITION 1M. Smith and D. V. Griffiths University of Manchester, UK JOHN WILEY & SONS Chichester . New York . Brisbane . Toronto . Singaporecme Copyright © 1982, 1988 by John Wiley & Sons Ltd. All rights rewerved. j 'No part ofthis book may be reproduced by any means, or transmitted, or translated ino machine language without the written permission of the publisher Library of Congress Cataloging-in- Publication Date: ‘Smith, LM. Programming the nite element method. Includes bibliographies and indexes. 1. Finite slement method —Computer programs 2. Engineering—Data processing. 3. Soil mechanics—Data processing 1. Grits, D.V. M1. "Tite. TASHTFSS6¢° 1988 620/001'515353. 87-8123, ISBN 0 471 91552 1 (cloth) ISBN 0 471 91553 X (paper) British Library Cataloguing in Publication Date: Smith, LM. an Moffat) Programming the finite element method. —2hd ed. 1,” Finite cement method—Data processing I Title “ Grifiths, D.V. 515353 TAMI.ES ISBN 0 471 91552 1 (loth) ISBN 0 471 91553 X (paper) ‘Typeset at Thomson Press (India) Limited, New Deli Printed and tound in Great Britain by Anchor Brendon Ltd, Tiptree, Essex tet eee Contents Preface to Second Edition Chapter 1 Preliminaries: Computer Strategies 10 Introduction 1.1 Computer Hardware 5 : 0 12 Store Management goobuG 13 Vector and Array Processors. =... . 5 14 Software. . 15. Portability of ‘Subroutines’ or ‘Procedures? 16 Libraries : Chapter 2 Spatial Discretisation by Fosse Elemene 20 Introduction 21 Rod Element... . 22 Rod Inertia Matrix. | ss 23. The Eigecvale Equation: - ‘i 24 Beam Element... 2. 1. 25 Beam Inertia Matrix. co 5 26 Beam with an Axial Force | rs 2.7 Beam on Elastic Foundation : 28 General Remarks on the Discretisation Process. | | 29 Alternative Derivation of Element Stiffness : 2.10 Two-dimensional Elements: Plane Stress and Strain 2.11 Energy Approach. . . as 2.12 Plane Element Inertia Matrix ©) |) | |) 243 Axisymmetric Stress and Strain : 2.14 Three-dimensional Stress and Sain...) |. |. 2.15 Plate Bending Element... : 216 Summary of Element Equations for Solids... | 2.17 Flow of Fiuids: Navier-Stokes Equations 2.18 Simplified Flow Equations . er 2.19 Simplified Fluid Flow: Steady State oo220 Siraplified Fiuid Flow: Transient State. 5 41 Conclusions. ©... oa 1a. 221 Simplified Fluid Flow with Advecton | clot 42. References Lt fell lll le 222 Further Coed Equations: Bit Consolidation | : 223 Conclusions. : Selle - Chapter 5 Static Equilibrium of Linear Elastic Solids... . 143 cee : $0 Introduction... . Dl as Chapter 3 Propanming Finite Bent Compattions Progam 30 se Sue Aste Ung Tacomas” 30 Introduction... a gles : af ance eee Progam 51: Pa St nas Ung Scale 3.2. Numerical Integration for Quadrilaterals . oe 33 Local Coordinates for Triangular Elements Prognim 52 Pane Sta Asay Using Fite rode Se ee Program 53 : Plane Strain Analysis Using Four-node 35. Element Assembly. : Plane Strain Analy “as 36 Incorporation of Boundary Conditions . Program 54 : Plane Strain Analysis Using Eight-node 37, Progamming Usig Buldig Boks Plane Stain Analy tet 38 phe Bor Rowger ON Progra $3: une Slt Asia sg i ede 3.10 Solution of Equilibrium Equations. |. Quadrilaterals . . . 165 2 a ns 311 Bralotion of Eigomales aod Eigeaecions | Propam 6: Auymeti Sia Asai Using Four 312 Solution of First Order Time-dependent Problems | Perit eee 33 Solution of Coupled Navier-Stokes Prosiens Sanne Une Beka Cuamexmmet in 3.14 Solution of Coupled Transient Problems.» a Pro oe . 2S Solution of Sond Onder Timedepndet Problems | | | Fee cde nim Acaivit Using Fournoded 3.6 Conclusions pit Program 59 : Threedimensional Analysis Using Bight RAT References. | |. pili Lapeer an 1 Chapter 4 Stati Bquiviam of Smuctes Progam 510: Sate, 40" Iowroducton a Sa Re . a ae Program 40 : Equilibrium of Uniform Beams. | oo danaG. Program 41 : Equilibrium of Stepped Beams Incorporating | Prescribed Displacement Boundary Chapter 6 Material Nor-linearity Se 12 Conditions. 60 Introduction. 192 Program 42 : Beam on Elastic Foundation with Numerically 641 Stress-Strain Behaviour | | | Poll l ll ase Integrated Siliness and Mass Matrices 62 Stessinvariants . | 195 Program 43 : Equilibrium of Two-dimensional Rigi oited 63 Failure Criteria Poll lll lasr Frames : 64 Generation of Body-oads | lS) 99 Program 44 : Equitbrium of Taree-dimensional Rigid- 65 Visco-pasticty inner’) jointed Frames... . 66 Initial Stress. 201 Program 4S : Equilibrium of Two-dimensional Pin-jointed 67 Comers on the Failure and Potential Suraces | | | | | 202 Troses Program 60 : Plane Strain Bearing Capacity Anaiysis— Program 46 : Equilibrium of Three-dimensional Pin-jointed Viscosplastc Method Using the Von Mises Truss : Criterion. 203 47 : Plastic Analyst of Two-dimensional Frames Program 61 : Plane Strain Slope Stability Ansiysis—Visco- Program 48 : Stability Analysis of Two-dimensional Frames plac Method Using the Moh-Covlonb 49 : Equilibrium of Rectangular Plates in Bending Criterion .Program 62 : Plane Strain Passive Earth Pressure ‘Analysis—Initial Stress Method Using the Mohr=Coulomb Criterion Program 63 : Axisymmetric ‘Undrained” Analysis —Visco- plastic Method Using the Mohr-Coulomb Criterion... be Program 6.4 : Three-dimensional Elasto-plastic Analysis— Visco-plastic Method Using the Mohr Coulomb Criterion 68 References. . ‘Chapter 7 Steady State Flow... Program 70 : Solution of Laplace's Equation Over a Plane ‘Area Using Four-node Quadrilateals. Program 7.4 : Ansys of Pana Fre Surface Flow Using Fournode Quadleteras 71 Reference Chapter 8 Transient Problems: First Order (Uncoupled) 80 Introduction Se oe : Program 80 : Solution of the Conduction Equation over a Rectangular Area . big Program &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 83 : Solution of the Conduction Equation over a Rectangle wing on EBE Product Algoritm Comparison of Programs 8.0,8.2 and 83. . . Program 8.4 : Diffusion-Convection Problem over a Rectangle —Transformed Analysis . Program 85 : Diflusion-Convection Problem over a Rectangle—Untransformed Analysis 82. References . Se Chapter 9 ‘Coupled’ Problems 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 « Rectangle in Plane Strain (Eight-node/Four-node Elements) 92 Referens ee nn Chapter 10 Eigenvalue Problems 100 Introduct Program 100 : Eigenvalues and Eigenvectors o! a Two” dimensional Frame of Beam-Column Elements—Consistent Mass Program 10.1 : Eigenvalues of a Rectangular Solid in Plane ‘Strain—Four-node Quadrilaterals with Lumped Mass... Program 102 : Eigenvalues and Figenvectors of a Rectangular Solid in Plane Stran—Fight- node Quadrilaterals with Lumped Mass Program 10.3 : Eigenvalues and Eigenvectors of a Rectangular Solid in Plane Stress—Eight- node Quadrieterals with Consistent Mass . Program 104 : Eigenvalues and Bigenvectors of a Rectangular Solid in Plane Stran—Four- node Quadrilaterals with Consistent Mass 101 References boo : Chapter 11 Forced Vibrations. 110 Introduction . Program 11.0 : Forced Vibration of a Rectangular Solid in Plane Strain Using Eight-node Quadrilaterals, Lumped Mass, Modal Superposition Method Program 11.1 : Forced Vibration of a Rectangular Solid in lane Strain Using Eight-node Quadrilaterals—Lumped or Coasistent Mass, Implicit Integration by Theta Method Program 11.2 : Forced Vibration of a Rectangular Elastic Solid in Plane Strain Using Fight-node Quadrilaterals—Lumped or Coasistent Mass, Implicit Integration by Wilson Theta Method Program 11.3 : Forced Vibration of a Rectangular Solid in Plane Strain Using Eight-node Quadrilaterals—Lumped Mass, Complex Response Method. Program 11.4 : Forced Vibration of a Rectangular Elasto- plastic Solid in Plane Strain Using Fight-node Quadriateras—Lumped Mass, Expt Integration . Program 11.5 : Forced Vibration of an Blastic Solid in Plane Strain Using Four-node Quadrilaterals Lumped or Consistent Mass, Mixed ExpliciyImplicit Integration LLL References... : . 31s - 319 320 = 320 320 325 330 333, 337 . 345 ~ 350Chapter 12 Analysis of Piles and Pile Groups 120 Introduction... ce Program 120 : ¢-z Analysis of Axially Loaded Piles Using ‘Two-noded Line Elements... . Program 124° p-y Analysis of Laterally Loaded Piles Using ‘Two-noded Beam Elements... . . Program 12.2 : Analysis of Vertically Loaded Pile Groups ‘Using Two-noded Line Elements—t-2 Springs for the Soil and Mindlin’s Equation for the Interactions : : Progsam 123 : Pile Drivability by Wave Equation—Two- node Line Elements—Implicit Integration in Time . ce 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 351 351 359 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 bas been dropped. However, the ‘geomechanics flavour remains, particularly in Chapters 6, 7, £, 9 and 12. Although the programming philosophy remains intact, Chapter 1 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 applic- ations and the equations of fuid low including the Navier-Stokes equations. Inthe original Chapter 3, choice of storage strategy was rather limited, s0 more versatility bas been introduced, in the form of a ‘profile’ (skyline) option, ‘Chapter 4 has been extended to embrace elastic, elasto-plastic and elastic stability analyses of all commonly encountered skeletal structures, while in ‘Chapter 5 a 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 probiems. In the present edition, transient problems are remosed 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 frst order transient problems. “The oniginal Chapter 8, now Chapter 9, has been extended to deal with two typical ‘coupled’ problems —those arising from the (steady state) Navier~Stokes ‘equations an¢ 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 eigeavalue problems is demonstrated. In the new Chapter 11 second xi‘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 singe 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 tae finite element 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 the computer 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 ‘or computations via the finite element technique. At the heart of what will be described is not a ‘program’ ‘or aset of programs but rather a collection (library) of procedures or subroutines Which perform certain functions analogous to the standard functions (SIN, SQRT, ABS, etc) provided in permanent library frm in all useful scientific computer languages. Because ofthe matrix structure of finite element formula- tions, most ofthe building block routines are concerned with manipulation of matrices. ‘The aim ofthe present book sto teach the reader to write intelligible programs and to use them, Super-fficiency 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 which itis 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 12 Table LL ‘Approximate Wordlength Core storage ct 1985 Designation Example (Gis) (words) Virtual memory? (€UK) Mico TBM PC 16 12a No 2K Super micro ICL PERQ 3 ax, Ye 2K Misi VAX 11/780 2 aK Yes 100K Mainframe CDC 7600 os 100K No 000K ‘Soper mainframe CDC 10GF 6 46000K Yes op00K book. In practice, hardware will ange from ‘micro? computers for small analyses and teaching purposes to ‘super’ computers for very large (especially non-linear) analyses. Itis a powerful feature of the programming strategy proposed that the same software will ran on all machine ranges. The special features of vector and array processors are described in a separate section (1.3) ‘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 ofa true subroutine structure in the 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 ofa ‘mainframe’. Thus a job taking two minutes on themainframe 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 isthe user is, advised not to operate at the extremes of the hardware’s capability. 12 STORE MANAGEMENT In the programs in this book it will be assumed that sufcient 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 can be of siz, say, 10,000 by 200. Thus a computer would have to have a main store of x 106 words to hold this information, and while some such computers ‘exist they are sill comparatively rare. A much more typical core size ist ofthe order of 1 x 10° words. ‘Two distinct strategies are used by computer manufacturers (o get round the problem of limited random access core. In the frst (simplest and earliest) 3 technique different ‘levels’ 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- fast large core memory’. Further back-up storage is available on a big disc but of ‘course when the programmer decides to transfer information to and from disc 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 queve. The second strategy removes store management from the user'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 on the ICL ATLASin 1961.Itcan be seen from Table 1.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, et.) fully utilised. (Clearly itis necessary forthe system to be able to translate the vrtual address of variables into a real address in core. This translation usually involves a complicated bit-pattemn 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. Hewever, store ‘management can never be totally removed from the user’s control. I: 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 108 Wordsina random fashion the paging requests wll readily ensure that very litle execution ofthe program can take place. 1.3 VECTOR AND ARRAY PROCESSORS Tn ater chapters it will be shown that finite element computations consist 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 bein the form of central processing ‘eetor 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, atthe very least, software can be written to pre-process standard high level language code to see ifit 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 ofthe special purpose software which has been written to date, for example to process arrays with various sparsity patterns, will become redundant, However, such 2 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. 1.4 SOFTWARE. Since all computers have different hardware (instruction formats, vector capability, ete) 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, etc) is ‘translated into the machine order code by a program called a ‘compiler’. tis convenient to group the major high level languages into those which are FORTRAN-ike (for example PASCAL) and those which are ALGOL-like (or ‘example ADA), The major difference between these groups that thelatter allow dynamic management of storage at run-time whereas the former do not. This difference has immediate implications for store management strategybecause itis clear that unless a computer permits genuine mult-programming and e one-level store in which programs are dynamically expanding and contracting there is not muuch 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, 50 one Would expect ALGOL-like languases to become more and more attractive. ‘There are many texts describing and comparing computer languages, often from the point of view of computer scientists, eg. 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 ALGOLdike and FORTRAN-like languages as far as the engincer-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 element rogram. In ALGOL 60 the program might read as follows: begin integer J, m,n: statements which read or compute I, m, 1; begin array o[:1.1-mi} BCI:m ln}, ofd:h.:m)r integer i,k; real x ‘statements which initialise a, b; for i:=1 step 1 until | do for j:= 1 step 1 until n do begin x:=.0; for k:=1 step 1 until m do x4 afi] EK: end; cli, =x statements which do something with end ead 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 DO11=1,L DO1J=1,N X=0 DO2K=1,M 2 X=X+AGKBKD CL=x 1 CONTINUE Statements which do something with C END Note that the declarations of the simple integers and the reali this program are unnecessary. Comparison of the two codes will easily show that the syntaxes of the two languages are very similar, eg. compare the ALGOL for t:= 1 step 1 until 1 do ‘with the corresponding FORTRAN DO1I=1,L fines 1 and 2 of ‘or the assignment statements + ali kJ+b(K, fi6 and X=X+AGK)+BK,D. 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 useless if we ‘want to multiply together arrays with dimension greater than 20 Itis this lack of ‘dynamic’ storage which is FORTRAN’s major limitation. It cam partially be alleviated by various devices such as allocating all storage to a large vector and keeping track ofthe starting addresses of the various arrays held in the vector, but cone can never get away without a dimensioning statement of some kind within a FORTRAN program. ‘This has various bad effects, among which is a 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 ofthe machine. Since most large machines allow multi-programming this can be inefficient 1S 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 rogram. 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 ia 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 right be as follows: procedure matmult (2, 6,¢l,m,m): valve fm, m; array a, b,c: integer I,m begin integer i,j,k; real x: for i:= 1 step 1 until Ido for j:=1 step 1 until x do jn x:=0 for k= step I until m do xx of KOK: end th Suppose we now wish to multiply the I x m) array a by the (m x n) array b to sive the (I n) array ¢ as before, and then to multiply c by the (n x ) arnay d to yield the (Ix D array e. The ALGOL 60 program would be as follows: begin integer J, m, m: procedure marmult (a,b, 6.747); read of compute J, m, m; i Lm, BEL, 1:n, eft sh din), dfn, 10, [1:11 statements which initalise a, ; matmult(2,b,6, mm) smatmul (dye, statements which do something with e end In the above, the procedure mavmultis assumed to be available in a user ibrary in exactly the same way as LOG, COS, etc, are available in the permanent library. The program structure can be seen 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 from our 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, - DO11=1,L DOIJ=1,N xX=0 DO2K=1,M 2 X=X+A(,K)*B(K,D), COD =X 1 CONTINUE RETURN END ‘Now the dimensioning statement inline 2 of this subroutine is merely adevice 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 cculd as ‘well have written REAL A(20,*), B(20, #), C(20,*), XMore importantly, since in FORTRAN the local A, B, C ate mapped onto the arrays of the main program, it can be seen that our subroutine is useless, unless arrays processed by i 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 js 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 ofthe 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, +), BOB, +}, CIC, +), X INTEGER IA, 1B, IC, L, M,N, I, J, K DO1I=1,L DOIJ=1,N X=0 DO2K=1.M 2 X=X+AGK)BKJ) C3) =X 1 CONTINUE RETURN END We could now construct our general matrix multiplication program asfolows, 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 IA/20/, 1B/5}, 1C/20/, 1D/1/, 1E/20/ Statements which read or compute L, M, N Statements which initialise A, B CALL MATMUL (A,IA,B, B,C, ICL, M,N) CALL MATMUL (C,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 lc onlyin the need to dimension arrays in the main program and ‘to pass the extra column size parameters in the subroutine calls. Infact we will -only be concerned in our fnite element programs with adjusting the sizes ofa few ‘arge arrays at runctime, and so the inconvenience can really be restricted to ‘changing a line at the beginning of the main program. This is most easily cone by the PARAMETER statement in FORTRAN 77. ‘Using this facility the second and third lines of the preceding program can be replaced by: PARAMETER (IA = 20, JA = 5, IB = 5, JB = 15, IC = 20, JC=15, *ID = 15, 3D =20, IE =20, 3E=20) REAL A(IA,JA), BOB, JB), CIC JO, DAD, JD), ECE, JE) so that only the PARAMETER line need ever be changed. This is the method used in this book. 16 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 preseat 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 infinite element analysis, and second it uses 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 be considered as “black box’ routines, that is the user need not understand or see them (how many users know how STN(X) is computed”). Scme ofthese 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 aso. ‘The test programs which are built up from the library routines are described in ‘Chapters 4 to 12and also listed there. Interested readers cansubstitute 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 subroutinecalls which could be substituted for the FORTRAN subroutines if required. ‘The programs described in this book are all assumed to be processed ‘in core’, Naturally, one could construct out of core versions ofthe library routines, which process large arrays for use 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 cur programs will be seen to be a nested structure and we will use representations called ‘structure ccharts' (Lindsey, 1977) rather than flow charts to describe their actions.10 ‘The main features of these charts are: () The Block This will be used for the outermost evel of each structure chart. Within a block, the indicated actions are to be performed sequentially. (i) The Choice ‘This corresponds to the if...then...else if...then...end if kind of construct. ii) The Loop “This comesin various forms, but weshall usually be concerned with the‘DO" loop. In particular, the structure chart notation encourages the use only of those GOTO statements that ate safe in a structured program. Using the notation, our matrix multiplication program would be represented as follows: u Imtiaz variables and roysandb ‘ow x colurm products Dosomething withe ‘The nested nature of a typical program can be seen quite clearly. 18 OTHER LANGUAGES ‘New programming languages appear ftom time to time, usually at the instigation ‘ofcomputer scientists rather than engineers. For exampie, PASCAL is popular as a teaching language and hasa simple structure enabling compact compilers to be written. In PASCAL 3 our matrix mukply 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12 FOR K:=1 TO M DO X+A(LK]+BEKN COX END; END; VAR L, M, N: INTEGER; ‘A:(1.10,1.10] OF REAL; Bi{1.10,1.10] OF REAL; Cif1.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, theres no real concept ofa stack pointer and in the main program A, B and Cneed to be assigned fixed dimensions. The PASCAL language does possess the built-in procedures (futctions) HIGH and LOW which enable bounds of arrays to be inspected; eg. 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 llanguage in which thisis possible is ALGOL 68. It has an extremely flexible structure and poiats the way towards the probable shape of future languages. Our matrix multiply program becomes: (proc matmul = ([,] real a, b, ref [,] real c) void: (int I= 1 upb a, m= 1 apb b, n=2 upb; b oc int i, j,k: for { to Ido for j ton do loc real x:= 0; for k to m do x+ = ali) bLk fod: els fAsm x od 0 {oe int 1m, read (m,n): Joe (121, 1:m) real a, loc [1:m, 1:n] real b, loc [1+l, 1:n] real ¢; read ((a,8): ‘matmul (a, b,¢), print (c) ) ‘The procedure calling parameters have been reduced to three (compared with FORTRAN’: nine) and the language has many more useful features. For B example, operations between matrices can be defined so that be dispensed with and operation c=axb 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 bad best be written in order to assure maximum readership and portability. (Theze are inevitably minor differences between machines, for example in their handling of free format input which is used exchisively 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, 2 library of subroutines san 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 vritten 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 forall the common machine ranges. ‘The structure of the remainder of the book is as follows. Chapter 2shows how the diflrential equations governing the behaviour of solids and fluids are programs are written ia a reasonably ‘structured? style, and tapes or discs of brary and programs will be supplied to readers on request. Versions are at present availabie 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, parti inthe authors’ field of geomechanics. However, the methods and programs described are equally applicable in many other fields such as structural mechanics, water resources, electiomagnetics and so on. Chapter 4 leads off with static analysis of skeletal structures and plate problems. Chapter 5 deals with static analysis of Tinear solids, whilst Chapter 6 discusses extensions to deal with material non- linearity. Chapter 7 is concerned with problems of uid low 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 sold and ffuid phases is treated, with applications to consolidation processes. A second type of ‘couplizg’ involves the Navier-Stokes equations. Chapter 10 contains Programs forthe 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 deseribed 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.10 REFERENCES Barron, D. W.(1977) An Introduction o the Study of Programming Languages. Cambridge University Press. ' Dijkstra, EW (1976) 4 Discipline of Programming. Prentice-Hall, Englewood Cis, New Jersey, Ford, B, and Sayers, D. K. (1976) Developing a single numerical algorithms library for diferent machine ranges. ACM Trans. Math. Software, 2, 115. Lindsey, C. (1977) Structure chart: 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 clements), This results io 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 deseribed in many texts, for example Zienkiewicz (1977), Strang and Fix (1973), Cook (1981), Connor and Brebbia (1976) and Rao (1982), but the principles will brielly 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 elastic rod, with end odes 1 and 2. The element has length L while u denotes the longitudinal displacements of points on the rod which is subjected to axial loading only. IFPis the axial force in the rod ata particular section and F is an applied body force with units of force per unit length then Pe od~BAe= BAM en and for equilibrium from Figure 2.100), ap Zireo 0 Hence the diferentil equation to be soived is ou eASS + Fao 2 1516 placement _prooy Force tb) te Figure 2.1. Equilibrium of a rod element (Although this isan ordinary differential equation, future equations, for example (2.3), preserve the same notation) In the finite element technique, the continuous variable «is approximated in terms ofits nodal values, u, and u,, through simple functions ofthe space variable called ‘shape functions’. That is ws Nyuy + Nats or vaty, waftbene 24 In simple examples the approximate equality in (24) could be made exact by oting x Bi m-(1-}), MX es) if the true variation of u is 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 e uy eZ, wai{t}+rao 26 ‘so that the partial differential equation has been replaced by an equation in the dliscretised space variables u, and u,. The problem now reduces to one of finding “good! values for uy and ws . ‘Many methods could be used to achieve this, and Crandall (1956) discusses ” collocation, subdomain, Galerkin and least squares techniques. Of these, Galerkin’s method, eg. Finlayson (1972) is the most widely used infinite element ‘work. The method consists in multiplying equation (2.6) by the shape functionsin turn and integrating the resulting residual over the element. Thus 2 IN, uy 2 ON, [fri}esson waftene [fii}eaeno on Note that inthe present example, double differentiation of the 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 Gntegration by parts) to yield typically [vZBlox= ~ [Fees + boundary terms which we ignore (28) ‘Hence (2.7) becomes aN, aN, ne Hi) [PEM aman, amram | {ant -FL tee Ox dx Ox x eye ie {ale ~{h 210) ‘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 aot EWE fu)_ fry cay feh-{ah aun) Tt which are the familiar ‘stiffness’ equations for a uniform prismatic rod, In matrix notation we may write RIE RIE KM ‘where KM is called the ‘clement stiffness matrix’, (212) 22 ROD INERTIA MATRIX Consider now the case of an unrestrained rod in longitudinal motion, Figure 2.1(¢) shows the equilibrium of a segment in which the body force is now8 ‘given by Newton's law as mass times acceleration. Ifthe mass per unit volume is », the differential equation becomes @u_ eu EAE AG = 0 2.13) 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, e My -[ffteadiew, waft bax ew « *TN,N, NiNz # fu -o4 (Tian man af} en Evaluation of integrals yields 4 aye fu -oaft JE} 19 7 o 4 z * Toe) . clement. (6) Plane general quadrilateral Figure 31 (@) Plane tangas clement.) oa > 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: ‘s0 that side 1 2 has ¢ = — 1, side 3.4 has { = +1, side 1 4 has y= — 1 and side 23 has y= + 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,), (1,1), (1, —1) respectively. . ‘When this choice is adopted, the shape functions for a four-noded quadrilateral with comer nodes take the simple form: Ny -H- 9-2) Na=KI- OU +9) em Ny=H+ Od +0) Nn 4 +90 =7) ‘and these can be used to describe the variation of unknowns such as displacement ‘or fluid potential ix an element as before. 9 ‘Under special circumstances the same shape functions can also be used to specify the relation between the global (x,y) and local (én) coordinate systems. If this is so the element is of a type called ‘isoparametric’ (Ergatoudis et al, 1968; Zienkiewicz et al, 1969), the four-node quadrilateral is an example. The ‘coordinate transformation is therefore Nyy + Naka + Nats + Nake Nx Naya + Nava + Naya + Nave =Ny 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, bestexpressed in matrix form by G2) (a) fox af) [2] ae EB) Jax Jax af fas allat fa 63) on} an an) (ay, oy, or ie} ay art _ 54 Ja a oJ jo (3.4) ay, 7 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: Jferor= f° * tial dgdn 6s 2 2 2 4 " , 7 a er Figure 33. (a) Degenerate quadrilateral. (6) Unacceptable ‘quadrilateralTable 3.1 Coordinates and weights in Gaus- sian quadrature formulae n cy ~My Hl ° 2 1 2 at 1 B v3 5 3 2 5 0 eI 5 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 (,) ‘become curvilinear. In practive (3.5) ate evaluated numerically, using Gaussian quadrature over quadrilateral regions (rons, 1966a, 19666). The quadrature rile in two dimensions are all of the form fi fi semacana $$ vonsteond 66) ‘where w, and w, are weighting coefficients and &,,n) are 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, ¢.g. Kopal (1961). The table assumes that the range of integration is +: 1; hence the reason for normalising the local coordinate system in this way. The approximate equality in (3.6) is exact for cubic functions when n= 2 and for quintics when n = 3, Usually one attempts to perform integrations over finite foo 4 0 t + As with quadrilaterals, numerical integration can be exact for certain polynomials. Forexample, in Table 32, the one-point rules exact for integration of frst 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 given by matrix equations, for example the equilibrium equation KP6=Q G12) 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 is to assemble the clements and so derive the properties ‘ofthe three-element system, Each element possesses node numbers, shown in ‘parentheses, which follow the scheme in Figures 3.1 to 33, namely rumbering ‘clockwise from the lower left hand corner. Since there is only one urknown at Figure35 Mesh of quadrilateral clemeats 3 ‘every node, the fluid potential, each individual element equation can be written KP, KP,s KPia (6) [Qs KP,2 KPzs KPza| }d:\ _ JQ: KPs2 KPs3 KPsal )osf~ \Q3 KP, KPss KPeal [s} [Oe However, in the mesh numbering system, not in parentheses, mesh node 4 corresponds to element node(1)in element 1 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 1 and term KP3, from element 3 would be added together ‘and would appear in location 4,4 and so on, The total system matrix for Figure 35 is given in Table 3.3, where the superscripts refer to element numbers. ‘This total system matrix is symmetrical provided its constituent matrices are symmetrical. The matzix also possesses the useful property of "bandednes? ‘which means that the terms are concentrated around the leading 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’ 1W 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 1’ = 13 compared with the other scheme with 1 = 2 If system symmetry exist 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 a maximum of four terms to the left or right ofthe diagonal in each row). Often, with a slight decrease of efficiency, the symmetrical half of a band matrix is stored as a rectangular array with a size equal to the G13 ‘Table 33 System stiffness matrix for mesh in Figure 35. Superscrips indicate element ‘mumbers KP; KPly 0 KP: KP, 00 KP KPh,+KPl, KPL, KP$, KPLa+KP3, KPLe 00 0 KP, KP, 0 KP3, KP 00 x KPL, 0.—OKPL,+KP2, KPL,+KP2, 0 KP2, KP. KPL, KPis+KPla KPLs KPL\+KPly KPLa+KPiy KPLe KPLi KPLs +KPLy o KPI, KPa ° KL, KPlg ° ° 0 KPa KPI ° ° ° 0 KP °Figure 3.6 Alternative mesh ‘numbering schemes rnumbe: 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 ‘on 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 freely floating elements or meshes are sometimes required but normally in eigenvalue 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 s singular ‘and this set of equations has no solution, ‘The simplest type of boundary condition-occurs when the dependent variable in the colution 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 to 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 102, to the leading diagonal ofthe stiffness’ matrix in the row in which the prescribed value is required. The term in the same row ofthe right hhand side vector is then set to the prescribed value multiplied by the augmented ‘sifiness' coeficient. For example, suppose the value ofthe uid head at node Sin Figure 35 is knowa to be 570units. 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 (Ky.s + 10"2)6, + small terms = 570 x (Ks, + 10) 14) which would have the eflot of making dy caual to $70, Cleaiy this procegures ‘only successful if indeed ‘small terms’ are small relative to 107, Boundary conditions can also involve gradients of the unknown in the forms G15) G16) ab ro 17) ‘where n is the normal to the boundary and Cy, C; are constants, ‘To be specific, consider a solution of the diffusion-advection equation (2.121) subject to boundary conditions (3.15), (3.16) and (3.17) respectively. When the second order terms c,(0*g/0x") and c,(@9/6y*) are integrated by parts, ‘boundary integrals of the type fort arise, where sis a length of boundary and |, the direction cosine of the normal. ary theca yan = O present no ditiaty soe the contour integral (3.6) vanishes. a Omer (16 gies roto an extra ntegral which forthe boundary lement jown in Figure 3. 6.18) f ENTC idl ds G19) ‘When ¢ is expanded as N@ we get an additional matrix 00 0 OJ ~Gola-xp] 0 2 10 229) 6 0120 620, 0000 ‘which must be added to the left hand side of the element equations.36 Figure 37 Boundary conditions involv- ing non-zero gradients of the unknown For boundary condition (3.17) the additional term is f eNTCI,ds 21 which is just a vector 22) 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 ofthe type $ = 0 or 6/n = O are the most ‘commen and are easily handled infinite element analyses. Thecases # = constant or a$/én = constant x ¢ are somewhat more complicated but can be appro- priately treated. Examples ofthe use ofthese types of boundary specification are included in the examples chapters. 37 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 props ation calculations Tt is anticipated that most users will eect to pre-compile all of the building blocks and to hold theae 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 hes been mede 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 are the basis of finite element computations. These special purpose routines are listed in Appendix 5. The black box routines should be ‘thought of as 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 2x 2) MVMULT (multiplies MATADD (adés) ‘TREEX3 (inverts 3 x 3) matrix by vector) MATINY (inverts N x N) MSMULT ‘enultiplies NULL3 (3d array) matrix by scalar) Similar routines associated with vectors are VECCOP VVMULT jcross-product) YECADD 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 spit into reduction and forward/backward re-substitution phases. BANRED CHOLIN GAUSBA BACKI — CHOBKI —SOLVBA BACK2 — CHOBK2 Method Gauss Choleski Gauss Equation Symmetrical Symmetrical Unsymmetrical ‘coefficients COMRED — SPARIN COMBKI —SPARBI COMBK2 —SPARB2 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 tigenvector determin- ation, for example58 Tridiagonalisation TRIDIA BANDRD Finds eigenvalues EVECTS BISECT Method QR factorisation Jacobi Equation coefficients Symmetrical Symmetrical, banded, fall Nx upper triangle stored stored ‘The first of each pair tridiagonalises the matrix and the second extracts all ofthe cigenvalues. The eigeavectors are found only by EVECTS. It should 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 case itis unlikely thatthe full range of eigenmodes ‘would be required. The various vector 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 LANCZI LANCZ2 are used to calculate the Lanczos vectors and eigenvectors of a matrix. Although simple matrix-by-vector multiplications can be accomplished by MYMULT, 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 LBBT whose action is described in Appendix 4. In a 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-dependent. Therefore siraple free format input to “TAPE 5 is adopted in all programs. The results of calculations are usually vectors or arrays and two subroutines enable the output ofthese in a standard E format to ‘TAPE 6. They are called PRINTV PRINTA respectively. 59. In order to describe the action of the remaining special purpose subroutines, which ar listed in full in Appendix 5, itis necessary first to consider the properties. of individual finite elements and then the representation of continua from assemblages ofthese elements. Static linear problems (including eigenproblems) are considered first. Thereafter modifications to programs to incorporate time dependence are added. 39 SPECIAL PURPOSE ROUTINES The job of these routines is to compute the element 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 ealeulation In the remairder of this chapter the notation adopted is that used in the subroutine listings in Appendix S. 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 clement stiffness matrix for plane elasticity given by (2.63) xan [forppacay 623) In program terminology tis becomes KM= f [ore'-pev-Beeaxcy 624 and ts formation is described by the inner loop of the structure chart in Figure 38. It is assumed for the moment that the element nodal coordinates (x;y) have been calculated and stored in the array COORD. For example, we would havefor « four-node quadrilateral Mao cooRD= 2 3 (3.25) Xe ‘The shape functions N are held in array FUN, in terms of local coordinates, asFindthe slementgeometr nodal coordinstes, COORD) ‘Nullthestitiness matrixspece, KM For allthe Gaussianintegrating points, NGF Find the Gauss pointcoordinates and weighting factorsin array SAMP. Formthe element shape functions, FUN and Formthe strain displacement matrix, BEE. Mutiply bythe stress strain matrix, DEE ‘Addthiscontribution tothe element ‘ifiness KM. ‘Assemble the current KM intothe global system matrix. Figure 38 Structure chart of element matrix assembly specified in (3.1) by 40 -XD(-ETA)|* 40 —XD0. + ETA) “)aq+XxDq+ETA) a 40 +XDG-ETA) ‘The B matrix contains derivatives ofthe shape functions and these are easily computed in the local coordinate system as oFUNF DER= |. 3 | or on FUN 6 (1-BTA) -(1+ETA) (1+ETA) ea ea (=X) (=X) (4X) 04x) a ‘The information in (3.26) and (3.27) for a four-node quadrilateral is formed by the subroutine FORMLN for the specific Gaussian integration points (XI, ETA), held in the array SAMP- where I and J run from 1 to NGP, the number of Gauis points specified in each direction. Figure 39(a) shows the typical layout and ordering for two-point Gaussian integration. In all cases SAMP is formed by the subroutine Gauss where NGP can take any value from 1 to 7. ‘The derivatives DER must then be converted into their counterparts in the (4,3) coordinate system, DERIV, by means of the Jacobian matrix transform- ation (3.3) or 3.4). From the isoparametric property, £} =COORD'+FUNT (628) and since the Jacobian matrix is given by ‘ax ay Re ax ay ann, itis clear that (3.29) can be obtained by differentiating (3.28) with respect to the local coordinates. In this way JAC = DER*COORD” 6.30) In order to compute DERIV we must invert JAC to give JACI using TWOBY?2 6.29) 4 b te |sah+ 7 1s ; Has + | Nopee \e 7 ‘aI (sy Figure 39 Integration schemes for (a) quadrilaterals and (b) trianglesa in this two-dimensional case and finally carry out the multiplication DERIV = JACI +DER 31) ‘Thus the sequence of operations CALL FORMLN (DER, IDER, FUN, SAMP, ISAMP, I, 3) CALL MATMUL (DER, IDER, COORD, ICOORD, JAG, AC, IT, NOD, IT) CALL TWOBY? (JAC, JAC, JACI, ACI, DET) 32) ‘CALL MATMUL (ACI, JACI,DER, IDER, DERIV, IDERIV, IT, IT, NOD) where NOD is the number of element nodes (four in tis 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 (%,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 Ctapter 1, IDER is the working size of array DER and so on, ‘The matrix BEE in (3.24) can now be assembled as it consists ofcomponents 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 four-node quadrilateral ELD = {Uy 0, Ma 0 thy 05 Uy D6} (3.34) ‘The variables u and vare simply the nodal displacements n the x and y directions respectively assuming the nodal ordering of Figure 3.5. ‘The components of the integral of BEE*DEE* BEE, evaluated at the Gauss points given by al combinations of | and J, can now be computed by transposing BEE to give BT, by forming the stress-strain matrix using the subroutine FMDEPS (for plane strain) hare 1 ° El-V) DEE-GSyaemliey } G35) 1eav 0 9 a= Vy, 68 and by carrying out the multiplications BIDB=BT+DEE*BEE 636 A plane stress analysis would be obtained by simply replacing FMDEPS by FMDSIG. ‘The integral is evaluated numerically by NOP Nor 2 2, WeWeBTDB,, 637) where W, and W, are the Gaussian weighting coefficients beld in array SAMP. ‘As soon as the element matrix is formed from (3.37) itis assembled into the elobal system matrix (or matrices) by special subroutines described later in this ‘chapter. Since the assembly process is common to all problems, modifications to the element matrix calculation for different situations will first be described, KM=DET+ 39.12 Plane strain (stress) analysts 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 FMTRIG ‘This delivers the shape functions (QL,—1L, 7 4L,L, QL, -DLy FUN= 4LLs (3.38) @L- DLs ails and their derivatives with respect to L, and Ly OFUNT ‘aL, DER=| ocunr OL, _f@u-) 4, 0 =, -@Ly-0) 45-1,) 0 Ly G@Ly~1) ALy-L,) -@ly-) 41, 639) 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 (Ly,L,) are held in the array SAMP and the corresponding weighting ‘coefficients in the vector WT. Both of these items are provided by the subroutine NUMINT ‘The version of this subroutine described in Appendix Sallows the total number ‘of integrating points (NIP) to take the value 1, 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, 1) CALL MATMUL (DER, IDER, COORD, ICOORD, JAC, UAC, IT, NOD, IT) G40) CALL TWOBY? (JAC, JAC, JACL, JACI, DET) CALL MATMUL, (GACI, TAC, DER, IDER, DERIV, IDERIV, IT, IT, NOD) (where 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=05eDET+ 5 eBTDB, Gan the factor of 0.5 being requited because the weights add up to unity whereas the area ofthe triangle in local coordinates only equals one hall A higher érder triangular element with fifteen nodes is also considered in ‘Chapter 5. The shape function and derivatives for this element are provided by routine FMTR1S5; otherwise the sequence of operations is virtually identical to 39.13 Axisymmmetric 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 (,), The stress-strain matrix is still given by an expression similar to (3.35) but the 4 x 4 DEE matrix is formed by FMDRAD In this case the integrated element stiffness is (2.66), namely ko= { [prepee+neeerdrde 049) 65 Considering a four-node element, the isoparametric property gives FUN”? = 1S Fun)-cooRD(K,) =SUM 643) Hence we have Napnae KM=SUMsDET+ W,+W,*BTD3,, G44) By comparison wit (3.37)it may be seen that when evaluated numerically the algorithms for axisymmetric and plane stiffaess formation will be essentially the same, despite the fact that they are algebraically quite diflerent. This is very significant from the points of view of programming effort and of program Aexibilty. However (3.44) now involves numerical evaluation of integrals involving 1/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 (i.e. SUM) approaches zero, Provided integration points do not lie on the r = Oaxis, however, reasonable results are usually achieved using a similar order of quadrature to planar analysis. 39.14 Plane steady laminar fluid flow 1 was shown in (2.115) that a fluid element has a ‘stiffness’ defined by xp [[reresay (45 ‘which becomes in program terminology KP= f f DERIV'+KAY+DERIV dx dy 646) and the similarity to 2.24) is obvious. The matrix DERIV simply contains the derivatives of the element shape functions with respect to (x,y) which were previously needed inthe analysis of solids and are formed by the sequence (3.32) ‘while KAY contains the permeability properties of the element in the form wave[E 8] om te compton, DTKD = DERIV+KAY*DERIV 3.48)66 is formed by appropriate matrix maltiplications and the final matrix summation for a quadrilateral element is KP =DET+'y) ). WieWy+DTKD, G49) 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,¢, (2.69), to take the general form Map { [renacay 850) where N are just the shape functions. In the case of plane fluid flow, since there is ‘only one degree of freedom per node, the ‘mass’ matrix is particularly simple in program terminology, namely Ma = RHO«[ [FUN sFUNA<¢y 631) By defining the product of the shape functions FTF = FUN'#FUN 652) we have nance =nnoete'S'S! WeW,eFTF, 53) MM =RHO+DETs 2 & 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 var=nwoe{[nemasey =RHO*DETs'S" 'S! WeW)2ECM,, 654) ‘When ‘lumped’ mass approximations are used MM becomes a diagonal matrix. For a four-noded quadrilateral (NOD = 4), for example, MM =(RHO*AREA/NOD)+I G55) where AREA isthe element area and I the unit matrix. For higher order elements, however, all nodes may not receive equal weighting. 67 39.16 Higher order quadrilateral elements To emphasise the ease with which clement types can be interchanged in, programs, consider the next member of the isoparametric quadrilateral group, namely the ‘quadratic’ quadrilateral with midside nodes shown in Figure 3.10. ‘The same local coordinate system is retained and the coordinate matrix becomes 1 cer xa coorp=|% 56 Xe Ye Or Xe Ya ‘The shape functions are now 41 —XN(1- ETA\(—X1- ETA -1) 40 -XD(1- ETA?) 4 —XD(0 + ETA)(—X1+ETA~1) 4U-XP)0+ETA) 4(1 + X(1 + ETAXI+ETA — 1) H1+XD(- ETA?) 4(1 + XI) — ETA)(XI- ETA -1) 40.- XI?) — BTA) ‘Which, together with their derivatives with respect to local coordinates, DER, are FUN- G57 Figure 3.10 General quadratic quadrilateral ‘element6 formed by the subroutine FMQUAD The sequence of operations: CALL FMQUAD (DER, IDER, FUN, SAMP, ISAMP, I,J) CALL MATMUL (DER, IDER, COORD, ICOORD, JAC, DAC, IT, NOD, IT) CALL TWOBY? (JAC, IAC, JACI, JACI, DET) G58) CALL MATMUL (ACI, DACI, DER, IDER, DERIV, IDERIV, IT, IT, NOD) (where NOD now equals 8) places the required derivatives with respect to (x, ))in DERIV and finds the Jacobian determinant DET. ‘Another plane element 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 4OXI(XI — 1)(BTA)(ETA — 1) 1 ~ HXD)QKL~ (ETA + (ETA ~ 1) HIT — s)(ETAMETA +1) 40K + 1X1 — N(ETAETA +1) FUN={ dexN(XT+ 1)ETA)ETA +1) 59) — HADI + ETA + META 1) AKDT + IMETAVETA 1) = H(XI + 1X1 = I(ETAMETA - 1) (X14 1)(K1 = 1)ETA + 1ETA~1) Figure 3.11 The Lagrangian nine-node clement Co) 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 different clement 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, FMLAG®, FMTRIG, ete, could be merged into a single routine). 39.11 Three-dimensional cubic elements ‘As was the case with changes of plane clement types, changes of element dimensions are readily made, For example, the eight-node brick element in Figure 3.12 i the three-dimensional extension of the four-noded quadrilateral Using the local coordinate system (7,0, the coordinate matrix is a Xt ane coorp=| ** %* (60a) Xe Ys 3s Xe Ye Xn Xe Ye te Figure 3.12 General linear brick element7” ‘The shape functions are 4(—X1(1 - ETA)(1 - ZETA) |* 4. —XD(1 - ETA)(I +ZETA) AU-+X00 — ETA) + ZETA) 4 +XD(1 — ETA)( - ZETA) $0. —XD(1 + ETA)(I - ZETA) 30 —XD(i + ETA)(L + ZETA) A. +XD( + ETA)(1 + ZETA) 4(.+X(1 + ETA) —ZETA) which together with their derivatives with respect to local coordinates are formed by the subroutine (3.606) FMLIN3 The sequence of operations: ‘CALL FMLIN3 (DER, IDER, FUN, SAMP,ISAMP, 1,J,K) ‘CALL MATMUL (DER, IDER, COORD, ICOORD, JAC, JAC, IT, NOD, IT) CALL TREEX3 (JAC, JAC, JACI, JACI, DET) G61) ‘CALL MATMUL GAC1, JACI, DER, IDER, DERIV, IDERIV, IT, IT, NOD) (where NOD now equals 8 and IT now equals 3) results in DERIV, the required sradients with respect to (x, 2) and the Jacobian determinant DET. ‘A higher order brick element with 20 nodes (Figure 3.13) is also used in 0 ttn 3 18 2 19 fl cs a7 Figure 3.13 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 stifiness is given by xat~ [{farepeeneedsdy ar 02) where BEE and DEE must now be formed by the subroutines FORMB3 FORMD3 respectively. The final summation becomes PoP ROP nor KM =DETe YY, Wie Wie Wee BIDB 663) where W;, W;, Wx are as usual Gaussian multipliers obtained from SAMP. For three-dimensional steady laminar fluid flow the element ‘stiffness’ is sr [ffosnov-xavspenw snc -{f DTKDaxdyaz oe, where KAY is the principal axes permeability tensor PX 0 0 KAY=/|0 PY 0 (3.65) Oo 0 PZ Gaussian integration gives NoPNr Nor KP=DETe YY, Wt WeWyeDTED yx 6.65) which is similar to (3.63). ‘The mass matrix for potential flow is wn.-nionfffrow-runsecne =RH0+[ frraxayde 657 which is replaced by quadrature as Np Nar NP MM=RHO+DET+ |) sae *WieWerFTFyx (3.68)n In solid mechanics, a three-dimensional equivalent of ECMAT would be required. This development is left to the reader. 39.18 Three-dimensional tetrahedron elements ‘An alternative element for three-dimensional analysis is the tetrahedron, the simplest of which has four corner nodes. The local coordinate system makes use ‘of volume ccordinates as shown in Figure 3.14. For example, point P can be identified by four volume coordinates (Ly,L3,Ls,L4) where G69) Asis to beexpected, one of these coordinates is redundant due to the identity L+hthtblat 670) ‘The shape functions for the constant strain tetrahedron are L|" Ly FUN= Ly G71) Ly 3 ' < * 2 Volume of 1 tetrahedron = Vy 3 Figure 3.4 A foursnode tetrahedron element B and these, together with their derivatives with respect toL,,Lzand L, are formed by the cubroutine FMTETS ‘The sequence of operations: CALL FMTETS (DER, IDER, FUN, SAMP, ISAMP,) CALL MATMUL_ (DER, IDER, COORD, ICOORD, JAC, JAC, IT, NOD, IT) CALL TREEX3 (JAC, JAC, JACI, JACI, DET) ‘CALL MATMUL, ‘GACI, ACI, DER, IDER, DERIV, IDERIY, IT, IT, NOD) 72) (where NOD now equals 4 results in DERIV which is needed to form the element matrices and DET which is used in the numerical integrations. The integration is performed by we KM=4sDET+ WyeBIDB, 67 ‘where Wis the weighting coefficient corresponding to the particular integrating point. The factor $s required because the weighting coefficients provided by the subroutine NUMIN3 always add up to unity whereas the volume of the tetrahedron in local ‘coordinates equals one-sixth. The addition of mid-side nodes results in the ten-noded tetrahedron which represents the next member of family. This element could easily be implemented by replacing FMTET by an appropriate subroutine. Transient, coupled poro-elastic transient and elastic-plastic analyses ll involve manipulations of the few simple element property matrices described above. Before describing such applications, methods of assembling elements and of solving linear equilibrium and eigenvalue problems must frst be described. 39.2 Assembly of elements ‘The nine special purpose subroutines READNF FORMKV FORMKU FMBIGK FKDIAG FORMKB FORMTB FORMKC FSPARV are concerned with assembling the individual element matrices to form the system matrix that approximates the desired continuum. Allied to these there must be a specification of the geometrical details, in particular the nodal ‘coordinates of each element and the element's place in some overall node numbering scheme.