Home

FEAP - - A Finite Element Analysis Program

1. 5 6 Elements with History Variables are aan aa sz ea an R en 5 6 1 Assigning amount of storage for each element 5 6 2 Accessing history data for each element HE 15 15 16 17 18 20 21 23 24 25 CONTENTS li 5 7 Elements with Finite Rotation Parameters 50 5 7 1 Nodal rotation treatment UROTmm subprogram 51 5 7 2 Local nodal rotation treatment uc posa Feo he Week 53 5 8 Energy Computation ria sola Ak da de A dille a dn Se ng 53 5 9 A Non linear Theory for a Truss an wee a arke LL 54 6 Utility routines 63 6 1 Numerical quadrature routines en ems a die 63 6 1 1 One dimensional quadrature e a SE a Mh a WA 63 6 1 2 Two dimensional quadrate qa a de Bak 0003 64 6 1 3 Three dimensional quadrature 65 6 2 Shape function subprograms ais Laan dd a 66 6 2 1 Shape functions in one dimension 67 6 2 2 Shape functions in two dimensions 68 6 2 3 Shape functions in three dimensions 69 6 3 Eigenvalues for 3 x 3 matrix ha oa R 70 Plot routines sar Ane E Libera s xudaya ved oya oek 71 Dala Mest pl ts rss ee a A adele 71 6 4 2 Element data plots iis pula yaz ae ode a cheek 73 6 4 3 S ani alias 76 List of Figures 4 1 Sample UMATIT modules waer yk SOE a ak oe ten a ik 19 4 2 Sample UMATLn module Lala h
2. values as v1 v2 Users may also provide additional input within the UMATIn module using the routines PINPUT or TINPUT described in Sect 2 1 The values of user parameters must be saved in the array ud the third argument of UMATIn In the current version there are 50 words of double precision values available by default Additional values may be allocated by assigning a larger value on the control record first record after the FEAP title record Each material model is assigned a user material number to the return parameter umat This number must be a positive integer Finally the number of history parameters to be assigned to each computation quadrature point must be returned in the parameter n1 Currently the parameter n3 may be set but is not available to the user material model Thus all history variables must be retained in the n1 list Use of history variables is described later as part of the UMODEL module CHAPTER 4 USER FUNCTIONS 19 subroutine umatil type vv ud n1 n3 rr Eu a a ee a ee 1 c Purpose User material model interface c Inputs c type Name of constitutive model character variable c vv Parameters user parameters from command line c Outputs c ni Number history terms nh1 nh2 c n3 Number history terms nh3 c ud User material parameters msn a printen nn n n n implicit none include iofile h logical pcomp character type 15 integer ni n3 real 8
3. Current values of the strains are as mentioned above passed in the array eps 6 and the trace of the strain in the parameter theta Thus 0 e 11 22 33 In addition if thermal problems are being solved the current value for the temperature is passed as td All material parameters for the current model are passed in the arrays d and ud The array d contains parameters assigned by standard FEAP commands as described in Tables 5 5 to 5 8 and the array ud contains values as assigned in the user module UMATIn For constitutive equations with additional internal variables that evolve in time users must define entries for the h1 array The number of entries available in the array for each evaluation i e each quadrature point is nh The value for nh is defined by the parameter n1 in module UMATIn see Fig 4 1 Values from the previous time step are passed back to the module in the array hn which also contains nh entries Users should never modify entries in the hn array Finally the values of the element operation switch is passed as the parameter isw See Chapter 5 for operations performed during different values of isw Using the above information users must compute values for the stress and the associated tangent matrix These are returned to the element in the arrays sig 6 and dd 6 6 In addition updates for any of the history parameters must be assigned in the array CHAPTER 4 USER FUNCTION
4. Name Description 61 K Fourier thermal conductivity 62 Ko Fourier thermal conductivity 63 K3 Fourier thermal conductivity 64 c Fourier specific heat 65 w Angular velocity 66 Q Body heat 67 Heat constitution added indicator 68 Follower loading indicator 69 Rotational mass factor 70 Damping factor Tl gi Ground acceleration factor 72 9 Ground acceleration factor 73 93 Ground acceleration factor 74 py Ground acceleration proportional load number 75 pa Ground acceleration proportional load number 76 ps Ground acceleration proportional load number 77 ao Rayleigh damping mass ratio TS ay Rayleigh damping stiffness ratio 79 Plate Shell Rod shear activation flag 80 Method Type 1 81 Method Type 2 82 Truss Rod quadrature number 83 Axial loading value 84 Constitutive start indicator 85 Polar angle indicator 86 Polar angle coord_1 87 Polar angle coord_2 88 Polar angle coord_3 89 Constitution transient type 90 41 Plane stress recovery 91 432 Plane stress recovery 92 a3 Plane stress recovery Table 5 7 Material Parameters 36 CHAPTER 5 ADDING ELEMENTS Parameter Name Description 93 sref Shear center type 94 y Shear center coordinate 95 y2 Shear center coordinate 96 lref Reference vector type 97 ng Reference vector parameter 98 na Reference vector parameter 99 n3 Reference vector parameter 100 Cross section s
5. The parameters defined in the common block are CHAPTER 2 DATA INPUT AND OUTPUT 6 ior input file unit number if negative input from keyboard iow output file unit number ilg solution log file unit number If an error occurs during input from the keyboard FEAP returns a value of true for the function and a user may reinput the record if the implied loop shown above is used For inputs from a file the program will stop and an error message indicating the type of error occurring and the location in an input file is written to the output file The input routines return data in a real 8 array td If any td i is to be used as an integer or real 4 quantity it must be cast to the correct type That is the following operations should be used to properly cast the variable type real 4 t real 8 td 5 integer j logical errck pinput errck pinput td 5 nint td 1 1 Integer assignment float td 2 I Real 4 assignment J t PINPUT may be used to input up to 16 individual expressions on one input record each input record is however limited to 255 characters The routine TINPUT differs from PINPUT by permitting text data to also be input It is useful for writing user commands or to input data described by character arrays The routine is used as logical errck tinput integer nt nn character text 16 16 real 8 td 16 errck tinput text nt td nn The parameter nt specifies the number of text values t
6. in the section of ELMTnn where isw 1 This parameter is located in the common erotas which has the structure real 8 xln real 8 rots rvel racc thkl integer rotyp common erotas x1n 9 9 4 amp rots 3 9 2 rve1 3 9 2 racc 3 9 2 thk1 9 rotyp The other entries in the common are arrays that return values for each element to treat the rotation values and rates We shall return to their description after describing the treatment of the global nodal data for rotations 5 7 1 Nodal rotation treatment UROTmm subprogram The nodal rotation data is stored in the array x1g which is dimensioned as x1g 9 6 numnp For node ng the entries in x1g are stored as follows CHAPTER 5 ADDING ELEMENTS 52 Description Rotation matrix A at time t Rotation matrix A at time tn Rotation matrix An41 at time tn41 Rotation increment angle A Rotation rate w at time t Rotation acceleration a at time tn Rotation angle 0 Rotation acceleration a at time t 41 Rotation acceleration aa at time trig Rotation matrix Ag at time to Component XLG 1 9 1 ng XLG 1 9 2 ng XLG 1 9 3 ng XLG 1 3 4 ng XLG 4 6 4 ng XLG 7 9 4 ng XLG 1 3 5 ng XLG 4 6 5 ng XLG 7 9 5 ng XLG 1 9 6 ng 1 hl Z O0 olo os While storage is provided for the 3 x 3 rotation matrices the representation may also be specified in terms of quaternions for which only 4 components are necessary In this case the 9 entries may be
7. real 8 common eltran real 8 common eluser integer common hdata integer common iofile integer real 8 common prstrs integer common sdata real 8 integer common tdata integer o head o head 20 numnp numel nummat nen neq ipr numnp numel nummat nen neq ipr nstep niter naugm titer taugm tform nstep niter naugm titer taugm tform iaugm iform intvc iautl nstepa iaugm iform intvc iautl nstepa dm n ma mct iel nel pstyp dm n ma mct iel nel pstyp tt tt 200 bpr ctan psil bpr 3 ctan 3 psil ut ut 200 nhi nh2 nh3 ht1 ht2 ht3 nlm plm nhi nh2 nh3 ht1 ht2 ht3 nlm plm ior iow ilg ior iow ilg nph ner int dshit nph ner erav j int 3 jshft ndf ndm nen1 nst nneq ndl nnlm nadd ndf ndm nen1 nst nneq ndl nnlm nadd ttim dt c1 c2 c3 c4 c5 chi dtcr idyn0 ttim dt c1 c2 c3 c4 c5 chi dtcr idyn0 n Up common pointer np 400 up 200 real 8 integer common hr mr hr 41 mr 1 Figure 5 2 FEAP Element Common Blocks CHAPTER 5 ADDING ELEMENTS include bdata h include cdata h include counts h include eldata h include elplot h include eltran h include hdata h include iofile h include prstrs h include tdata h include pointer h include comblk h Figure 5 3 FEAP Element Common Blocks using Includes Figure 5 4 2 Node Truss Element CHAPTER 5 ADDING E
8. Consistent Mass DAMPn n 16 compro Damping JPn n 20 neq Profile pointer LMASn n 12 neq Lump Mass TANGn n maxpro Symmetric tangent UTANn n 4 maxpro Unsymmetric tangent Table 3 2 Solution Array Names Numbers and Sized NAME Num dim 1 dim 2 dim 3 Description ANGL 46 nen Angle LD 35 nst Assembly nos P 35 nst Element vector 36 nst nst Element matrix TL 39 nen Temperature UL 41 ndf nen 6 Solution array XL 44 ndm nen Coordinates Table 3 3 Element Array Names Numbers and Sizes 11 CHAPTER 3 ALLOCATING ARRAYS NAME Description IX 1 e Global node 1 SER to IX nen e Global node nen IX nen 1 e IX nen 2 e IX nen 3 e IX nen1 e IX nen1 1 e IX nen1 2 e IX nen1 3 e H1 history data pointer H2 history data pointer H3 history data pointer Element material type number Element region number default 0 Active deactive indicator Active deactive start Table 3 4 Element connection array IX use for element e NAME Description IE 1 ma Global DOF 1 to IE ndf ma Global DOF ndf IE nie ma IE nie 1 ma IE nie 2 ma IE nie 3 ma IE nie 4 ma IE nie 5 ma IE nie 6 ma IE nie 7 ma IE nie 8 ma Number history variables element NH1 and NH2 Element material type number ELMTO1 1 etc Element material type identifier default ma Off
9. FLG where Parameter Description 56 2 Natural coordinates n XL NDM NEL Element coordinates in local order NDM Spatial dimension mesh 2 or 3 NEL Number nodes on element 4 9 12 16 IX NEL Element global node numbers FLG Return 7 derivatives if true or x y derivatives if false SHP 3 NEL Shape functions and derivatives XJAC Jacobian transformation from x y to n CHAPTER 6 UTILITY ROUTINES 69 The array SHP stores the values in the order SHP 1 i derivative with respect to or x SHP 2 i derivative with respect to y or y SHP 3 i shape function Two dimensional Co isoparametric interpolation on triangles of linear quadratic and cubic order may be obtained using the subprogram call CALL TRISHP SS XL NDM IORD XJAC SHP where Parameter Description SS 3 Area coordinates Li Lo Ls XL NDM Element coordinates in local order NDM Spatial dimension mesh 2 or 3 IORD Order of interpolation 1 3 node 2 6 node 3 7 node 4 6 node 3 bubble 10 10 node cubic XJAC Jacobian transformation from x y to n SHP 3 NEL Shape functions and derivatives The array SHP stores the values in the order SHP 1 i derivative with respect to or x SHP 2 i derivative with respect to 7 or y SHP 3 i shape function The parameter IORD defines the order of interpolation If it is 1 simple 3 node triangles with linear interpolation is re
10. and develop an algorithm to solve a possibly non linear problem In FEAP the integration method for nodal quantities is taken as a one step algorithm with each quantity defined only at discrete times t Accordingly we have displace ments u t with velocities and accelerations denoted as u tn Z Vi tn and ta a tn A typical example for an integration algorithm for these discrete quantities is New mark s method where 1 CHAPTER 5 ADDING ELEMENTS Parameter Name Description 1 E Young s modulus Di Poisson ratio 3 a Thermal expansion coefficient 410 Mass density 5 Quadrature order for arrays 6 Quadrature order for outputs Tia Mass interpolation a 0 Diagonal a 1 Consistent 814 Loading intensity plates shells 9 To Stress free reference temperature 10 lk Shear factor plates shells beams 11 by Body force volume in 1 directions 12 bo Body force volume in 2 directions 13 bs Body force volume in 3 directions TA Thickness plates shells 15 nhi History variable counter 16 stype Two dimensional type 1 plane stress 2 plane strain 3 axisymmetric 17 etype Element formulation 1 displ 2 mixed 3 enhanced 18 dtype Deformation type lt finite gt small 19 tdof Thermal degree of freedom 20 imat Non linear elastic material type 21 di Material moduli 22 das Material moduli 23 da Material moduli 24 dio Material modul
11. and residual the unused positions need not be defined the tangent array s and the residual r are set to zero before each element routine is called The arrays 1 1 ul i j 1 ul i j 4 and ul i j 5 described in Table 5 1 are used to obtain the nodal coordinates displacements velocities and accelerations respectively When FEAP solves a problem without transient loading e g inertial CHAPTER 5 ADDING ELEMENTS Al Parameter Description ctan 1 cj Multiplier of s matrix for ul i j 1 terms e g stiffness matrix multiplier ctan 2 c Multiplier of s matrix for ul i j 4 terms e g damping matrix multiplier ctan 3 cz Multiplier of s matrix for ul i j 5 terms e g mass matrix multiplier Table 5 9 Tangent Parameters loading as mass times acceleration the velocities and accelerations are set to zero prior to calling the element subroutine Consequently in programming the steps to compute the residual r the inertia terms have no effect for static or quasi static prob lems and may be included generally there are very few additional operations involved to add these terms The programming of the tangent array however must distinguish between cases in which transient e g inertial loads are present and those in which they are omitted The different cases are implemented in FEAP by making appropri ate assignments to the c parameters To facilitate the programming of the tangent array returne
12. as the first field For example in a solid element the command sequence can be MATErial ma SOLId UCONstitutive xxxx v1 v2 The role of the xxxx and vi data will be described in Section 4 2 1 CHAPTER 4 USER FUNCTIONS 18 It is possible to use standard input parameters defined in Tables 5 5 to 5 8 as well as by preceding the UCON command with a normal input sequence For example if isotropic elastic properties are needed they may be included in the input sequence as MATErial ma SOLId ELAStic ISOTropic e nu UCONstitutive xxxx vi v2 No standard commands should follow the UCON command Alternatively users may input elastic properties as part of their UMATIn module If the user routine does input additional data records after the UCON record and these are terminated by a blank record a second blank record will be needed to discontinue material data input for this set In all cases at least one blank record is always needed to terminate the input of standard options for the material set 4 2 1 The UMATIn Module A sample module for a user constitutive model is shown in Fig 4 1 As shown in this figure the UMATIn module has 5 arguments The name of the constitutive equation to be described is passed in the first parameter type The second parameter passes an array vv which may be used to define up to 5 parameters for the material model The example shown above for the UCON includes the type data as xxxx and the array vv
13. ensure that a calling sequence is not sensitive to a change in pointer One way pointer changes can still lead to errors is through a program call subname hr np 111 mr np 112 and then change the length of the array number 111 or 112 in the subroutine Chapter 4 USER FUNCTIONS Users may add their own procedures to facilitate additional mesh input features to per form transformations or manipulations on mesh data to add new solution commands or to add new plot capabilities 4 1 Mesh Input Functions UMESHn To add a mesh input command a subprogram with the name UMESHn where n has a value between 0 and 9 must be written compiled and linked with the program The basic structure of the routine UMESH1 is subroutine umesh1 uprt c User defined routine to input mesh data to FEAP implicit none include umaci h 1 Contains UCT variable logical uprt Set name mesl to user defined if pcomp uct mesi 4 then uct Set user defined command name else C User execution function statements follow end if 15 CHAPTER 4 USER FUNCTIONS 16 end The parameter UPRT is a logical parameter which is set to false when the NOPRint mesh command is given and to true when the PRINt command is used default is true The common block UMAC1 transfers the character variable UCT for the name of the command The default names are MESn where n is the same as the routine name number
14. member element shown in Figure 5 4 are the balance of momentum equation O Aoss Os the strain displacement equation for small deformations Ab pA SS Os and a constitutive equation For example considering a linear elastic material the constitutive equation may be written as Oss E ss Boundary and initial conditions must also be specified to obtain a well posed problem however our emphasis here is the derivation of the element arrays associated with the above differential equations In the above e s is the coordinate along the truss member axis e b is a loading in direction s per unit length e A is the truss cross section area e pis the mass density per unit volume e u is a displacement in direction s e is an acceleration in direction s v e css is a strain along the truss member axis and e Gss is the stress on a truss cross section CHAPTER 5 ADDING ELEMENTS 43 The equations may also be deduced from the variational equation d d ou Fest Ads S suipai ds Y f busts lLa L i 1 YL YL where dII contains the boundary and loading terms not associated with an element Where in addition to previously defined quantities we define e dis the spatial dimension of the truss 1 2 or 3 e x are the Cartesian coordinates in the d directions e L is the length of the truss member e du is a virtual displacement in direction x e is an accel
15. s i i s i j s j i s j j cmd cmo cmo cmd s i i s i j s j i s j j end do endif Figure 5 7 Truss Tangent Matrix Part 2 CHAPTER 5 ADDING ELEMENTS 62 subroutine slcn2d sig shp xsj sg lint nel nes p s gt a TAR eee c Purpose Project element variables to nodes c Inputs sig nes Stresses at quadrature points c shp nel Shape functions at quadrature points c xsj Volume element at quadrature points G sg 3 Gauss points 1 2 and weights 3 c lint Number of quadrature points c nel Number nodes on element c nes Dimension of stress array c Outputs p nen Weights for lumped projection c s nen Integral of variables es be 000 5 o o 1 implicit none include cdata h 1 Contains nen include strnum h 1 Contains iste integer i l lint nel nes realx8 xg p s men xsj sig nes shp nel sg 3 9 do 1 1 lint do i 1 nel xg shp i 1 xsj 1 pG pi xg s i 1 s i 1 sig 1 1 xg s i 2 s i 2 sig 2 1 xg s i 3 s i 3 sig 3 1 xg s i 4 s i 4 sig 4 1 xg end do i end do i iste 4 Returns number projections end Figure 5 8 Element variable projection routine Chapter 6 UTILITY ROUTINES The FEAP system includes several subprograms that can assist developers in writing new modules In the next sections we describe some of the routines which perform numerical inte
16. the arrays using their pointers e g in a subroutine call as include pointer h include comblk h call subname hr np 111 mr np 112 Note the use of hr and mr for the double precision and integer references respectively Also the use of the pointers avoids a need to include the array reference until it is needed in a computation A short list of the mesh arrays available in FEAP is given in Table 3 1 for solution arrays in Table 3 2 and for element arrays in Table 3 3 The names of all active arrays in any analysis may be obtained using the SHOW DICTionary solution command Tables 3 4 and 3 5 describe the use of individual entries in the arrays IX and IE respectively See the subprogram palloc f in the program directory for the names and numbers of existing arrays CHAPTER 3 ALLOCATING ARRAYS NAME Num dim 1 dim 2 dim 3 Description ANG 45 numnp Angle D 25 ndd nummat Material parameters F 27 ndf numnp 2 Force and Displacement ID 31 ndf numnp 2 Equation nos IE 92 nie nummat Element control dofs etc IX 33 neni numel Element connections T 38 numnp Temperature U 40 ndf numnp 3 Solution array VEL 42 ndf numnp nt Solution rate array X 43 ndm numnp Coordinates Table 3 1 Mesh Array Names Numbers and Sizes NAME Num dim 1 dim 2 dim 3 Description CMASn n 8 compro
17. vv 5 ud c Specify type of user model if pcomp type mat1 4 then type E 1d Specify new name for model c Input output user data and save in ud array c Set values of ni if required ni write iow User Constitutive Inputs E vv 1 ud 1 vv 1 endif end Figure 4 1 Sample UMATI1 module CHAPTER 4 USER FUNCTIONS 20 4 2 2 The UMATLn Module A sample for the UMATL1 module is shown in Fig 4 2 This subprogram will be called by many of the elements included within FEAP if a user model has been specified as part of the MATE mesh data see previous subsection The user model will not be called for truss frame plate and shell elements which use resultant models to describe behavior Also any form which requires a one dimensional model will not use a UMATLn module The module is designed to compute three dimensional constitutive models in which the stress and strain are stored as 6 component vectors and the tangent moduli as a 6 x 6 matrix Strains are passed to UMATLn in the argument array eps 6 and stored in the order T e en 22 33 712 723 931 where 2 is the engineering shearing strain Stress and moduli are to be associated with the same ordering and returned in the argument arrays dimensioned as sig 6 and dd 6 6 respectively All values are to double precision i e REAL 8 When UMATLn is called the model n will be that which is defined in the module UMATIn
18. 16 O J integrals Soln JINT 17 O Set after activation Soln ACTI 18 O Set after deactivation Soln DEAC 19 NOT AVAILABLE used in modal base BASE Uses isw 5 in element 20 O Element plotting Plot PELM 23 O Compute element loads only ARCL 25 O Zienkiewicz Zhu projection Soln ZZHU 26 R Used to compute mesh boundary Called by default Table 5 2 Task Options for FEAP Element Subprogram R Required O Optional E For eigensolutions CHAPTER 5 ADDING ELEMENTS 31 select case isw case 1 c Input material parameters case default end select the CASE DEFAULT should not perform any operations unless all options are programmed Finally if the form go to 1 2 isw return is used the RETURN statement should always be included as shown This prevents any unexpected execution of a statement that appears after the GO TO Some of the options for additional data passed through common blocks is shown in Figure 5 2 with each variable defined in Table 5 3 Also in Figure 5 3 the reference to common blocks using include statements is shown In the prototype routine the number of nodes on an element nen which is used to dimension ul is passed in the labeled common cdata Additional discussion is given below on use of some of the other data passed through the common blocks 5 1 Material property storage The material parameters to be stored in the array D with pointer np 25 may be input using the subprogra
19. Assignment of a unique character name which must not conflict with names already assigned for mesh input commands should be used to replace the xxxx shown When FEAP begins execution it scans all of the UMESHn routines and replaces the command names mest etc by the user furnished names Thus when the command HELP is issued while in interactive MESH mode the user name will appear in the list instead of the default name note FEAP does not always display all available commands To see all commands issue the command MANUal 3 and then the HELP command The ability to get array names as shown in Chapter 3 can be used to develop user routines for input of coordinates element connections etc With this facility it is possible to develop an ability to directly input data prepared by other programs which may be in a format which is not compatible with the requirements of standard FEAP mesh commands 4 1 1 Command line data It is possible to include extra data on the command line for user functions If the data is not passed as an argument e g as in the UMESHn functions the data is recovered from a character string yyy 255 contained in the labeled common statement chdata The string is formatted into words with a field widths of 15 characters Text words are left justified and numerical values are right justified The data must be recovered in the function before any additional input statements are processed For example if a user input functi
20. FEAP A Finite Element Analysis Program Version 8 2 Programmer Manual Robert L Taylor Department of Civil and Environmental Engineering University of California at Berkeley Berkeley California 94720 1710 E Mail rlt ce berkeley edu March 2008 Contents Introduction 1 1 Setting Program Options gt a a oa 1 2 Uses of Common and Include Statements Data Input and Output 2 1 Parameters and Expressions 2 ne oi Bane ed Scere a ears oe 2 2 Array Q tp tg A NP GATA y s Diani Allocating Arrays User Functions 4 1 Mesh Input Functions UMESHn 4 1 1 Command 4 4 2 User Material Models 4 2 1 The UMATIn Module 4 2 2 The UMATLn Module 4 2 3 Auto time step controls bara rt 4 3 Mesh Manipulation Functions UMANIn 4 4 Solution Command Functions UMACRn 4 5 Plot Command Functions UPL Tn Adding Elements 5 1 Material property storage suo tra tl as i Geet 5 2 Non linear Transient Solution Forms 5 3 Example 2 Node Truss Element co 20 KAN KANG 5 3 1 Theory for a Truss La ha PNG 5 4 Additional Options in Elements 54 1 Task Options eta di a ok AE oh 54 2 Task 03 TO s 7 5 5 Projection of element variables to nodes
21. LEMENTS 59 if isw eq 3 or isw eq 6 then c Compute element length L2 0 0d0 do i 1 ndm L2 L2 1 4 2 xl i 1 2 end do L sqrt L2 c Compute strain displacement matrix bb i 1 x1Gi 2 x1G 1 Lr bb i 2 bb i 1 eps eps bb i 2 ul i 2 1 ulGi 1 1 end do c Compute mass terms cmd rhoA L 3 0d0 cmo cmd 0 5d0 Form body inertia force vector dm prop ld sigA EAxepsxL body 0 5d0 L dm do i 1 ndm r i 1 body d 6 i bb i 1 sigA amp emd ul i 1 5 cmo ul i 2 5 r i 2 body d 6 i bb i 2 sigA amp cmo ul i 1 5 cmd ul i 2 5 end do Figure 5 5 Element residual for two node truss CHAPTER 5 ADDING ELEMENTS if isw eq 3 then c Compute element length L2 0 0d0 do i 1 ndm L2 12 1 4 2 x1 i 1 xx2 end do L sqrt L2 c Form stiffness multiplier dd ctan 1 FA L c Compute strain displacement matrix Lr 1 d0 L2 do i 1 ndm bb i 1 x1 1 2 xl i 1 Lr bb i 2 bb i 1 db i 1 dd bb i 1 db i 2 db i 1 end do Figure 5 6 Truss Tangent Matrix Part 1 60 CHAPTER 5 ADDING ELEMENTS 61 Compute stiffness terms ndm 4 or ndf il O do ii 1 2 jl 0 do jj 1 2 do i 1 ndm do j 1 ndm s itil j j1 db i ii bb j jj end do end do ji ji ndf end do 11 il ndf end do c Compute mass terms and correct for inertial effects cmd ctan 3 rhoA L 3 0d0 cmo cmd 0 5d0 do i 1 ndm j i ndf
22. S 21 hi and returned to the element Values of history variables returned are not used for all values of isw e g when reporting or projecting stresses under isw 4 and isw 8 they are not saved Values retained in the h1 array are copied to the hn array each time the command statement TIME is issued in a solution 4 2 3 Auto time step control The solution command AUTO MATErial rvalu 1 rvalu 2 rvalu 3 initiates an attempt to control the solution process by a variable time stepping algo rithm based on a user set value in the material constitution The value to be set is named rmeas which is passed between constitution and solution modules in the labeled common real 8 rmeas rvalu logical aratfl common elauto rmeas rvalu 3 aratfl The three parameters may be used in defining an acceptable value for rmeas The algorithm coded monitors the solution during a standard iteration process set by for example LOOP n TANG 1 NEXT If during any iteration up to n the value of rmeas exceeds a value of 2 rmeas 0 at the start of the loop a new value of At is immediately set to Atnew 0 85 At rmeas and the iteration process is started over On the other hand if convergence occurs during the time step and the value of rmeas is smaller than 1 25 the time step is adjusted according to Atnew 1 50 At rmeas lt 0 5 Aben 1 25 At 0 5 lt rmeas lt 0 8 At rmeas 0 8 lt rmeas CHAPTER 4 USER FUNC
23. TIONS 22 000000000000 aaa subroutine umatl1 eps theta td d ud hn hi nh istrt sig dd isw Purpose User Constitutive Model Input eps Current strains at point small deformation Deformation gradient at point finite deformation theta Trace of strain at point Determinant of deformation gradient td Temperature change d Program material parameters ud User material parameters hn nh History terms at point t_n hi nh History terms at point t_n 1 nh Number of history terms isw Solution option from element Output sig 6 Stresses at point dd 6 6 Current material tangent moduli 5 unu su gunu nunun n n n 1 implicit none integer umat nh isw i real 8 td real 8 eps theta d ud hn nh hi nh sig 6 dd 6 6 Dummy model sig ud 1 eps do i 1 6 dd i i ud 1 sig i ud 1 eps i end do endif end Figure 4 2 Sample UMATLn module CHAPTER 4 USER FUNCTIONS 23 Finally if convergence does not occur with in the n steps then the time step is reset according to Atnew 0 85 At rmeas 1 25 lt rmeas Atnew AY otherwise After any of the above adjustments the value of rmeas is reset to zero 0 An optimal value of rmeas is 1 25 which leaves the step unchanged The above algorithm was proposed by Weber et al 1 4 3 Mesh Manipulation Functions UMANIn The UMANIn modules where n ranges from 0 to 9 may be
24. alues are obtained using the call CALL PSTR3D SIG PV where SIG 6 stores stresses in the order 011 022 033 012 023 031 and returns principal values in PV 3 in the order 01 02 03 Roots are ordered from most positive to most negative 6 4 Plot routines Several options exist in the FEAP system to create graphical plots for data and results 6 4 1 Mesh plots FEAP has plot capabilities to represent some standard element shapes For continuum elements where the shape of the element is identical to the space dimension of the mesh i e NDM it is not necessary to provide any extra descriptions However if the dimension of the element topology is different from the mesh dimension it is necessary to add the include statement include eldata h CHAPTER 6 UTILITY ROUTINES 72 and the statement pstyp pdim within the ISW 1 part of the element routine Here pdim is the dimension of the element topology For example in a three dimensional shell problem where NDM 3 and the element topology is two dimensional the statement is given as pstyp 2 Provided the nodal numbering of an element is as described in the FEAP User manual i e numbered with vertex nodes first followed by mid side nodes then face nodes and finally internal nodes the program can use the actual number of nodes on the element to draw each element Failure to include a pstyp statement will usually result in unpredictable plots of the mesh and
25. always define the address of an array in a calling statement to also be 64 bits in length For example use of integer ioff ioff np 111 numnp call submat hr ioff would cause an error since the pointer ioff is only 32 bits in length To avoid this problem it is necessary to either declare ioff to be 64 bits long as integer 8 ioff or use one of the FEAP include files p int h defining the integer type array fp 10 or p point h defining the integer type scalar point A type of correct length is controlled at compile time by which subdirectory is used i e integer4 or integer8 Using this scheme permits direct reference to either real 8 or integer arrays in pro gram modules without need to pass arrays through arguments of subprograms A subprogram PALLOC controls the allocation of all standard arrays in FEAP defined by the np pointers and a subprogram UALLOC permits users to add allocation for their own arrays defined by the pointers up The basic use of the routines is provided by an instruction setvar palloc number NAME length precision or setvar ualloc number NAME length precision where setvar palloc and ualloc are logical types number is an integer number of the array NAME is a 5 character name of the array length is the number of words of storage needed for the array and precision is the type of array to allocate 1 for integer and 2 for real 8 types Upon initial assignment of any arr
26. antities it is necessary to add the pro jection operations for each element under ISW 8 These are performed locally for each element similar to all other operations Figure 5 8 shows a simple routine for two dimensional elements with 4 stress components begin projected When multiple element types are used in an analysis users must be careful to project like quantities to common values of the ST nen x array so as to get correct results Also when results are displayed it is necessary to plot results by material type to obtain correct indications of stress discontinuities at material interfaces 5 6 Elements with History Variables FEAP provides options for each element to manage variables which must be saved during the solution These are history variables and are separated into three groups a Variables associated with the last converged solution time t b Variables associated with the current solution time tn4 1 and variables which are not associated to any particular time All history variables are associated with the allocation name H which has a pointer value 49 Users are not permitted direct access to the data stored as H of course it is possible to access from hr np 49 but this should not normally 3An implementation of the Zienkiewicz Zhu projection method is implemented using ISW 24 CHAPTER 5 ADDING ELEMENTS 49 be attempted Before calling the element routine for each element FEAP transfers the required history v
27. arameters are identical to those for MPRINT except the array must be of type integer Chapter 3 ALLOCATING ARRAYS Dynamic data allocation is accomplished in FEAP by defining addresses in pointers contained in the common block defined in pointer h This common block contains pointers np for standard program arrays and up for user defined arrays and has the form integer num_nps hum_ups parameter num_nps 400 num_ups 200 integer np up common pointer np num_nps up num_ups Each pointer is an offset relative to the address of a REAL 8 array hr 1 or an INTEGER array mr 1 defined in a blank common real 8 hr integer mr common hr 1 mr 1 which is placed in the file comb1k h in the include integer4 or include integer8 directory The pointers in the integer4 subdirectory have 32 bit lengths and should be used for cases where addressing is within 4 GByte range The pointers in the integer8 subdirectory have 64 bit lengths and should be used for additional index ranges The arrays hr and mr are used to establish addresses only and not to physically store data This mechanism permits references to elements in arrays which have positions relative to hr or mr that may be after or before 1 Thus FEAP must be compiled without strict array bound checking Size of problems is limited only by the available memory in the computer used CHAPTER 3 ALLOCATING ARRAYS 9 When using 64 bit pointers users must be careful to
28. ariable to a local storage for each type Users may then access the history data for each element and if necessary update values and return them FEAP Only for specific actions will the local history data be transferred back to the appropriate H locations The element history data associated with t starts at the memory address of the pointer for NH1 using the double precision dummy array HR in blank common similarly data for tn41 starts at the memory address of the pointer for NH2 and that not associated with a time at NH3 The three pointers are passed to each element routine in the labeled common integer nh1 nh2 nh3 common hdata nh1 nh2 nh3 5 6 1 Assigning amount of storage for each element The specification for the amount of history information to be associated with each material set is controlled in the isw 1 task of an element routine For each material type specified within the element routine a value for the length of the NH1 and the NH3 data must be provided the amount of NH2 data will be the same as for NH1 This is accomplished by setting the variables nh1 and nh2 in common hdata see above to the required values That is the statements required are if isw eq 1 then nhi nh3 6 10 reserves 6 words of NH1 and NH2 data and 10 words of NH3 data for each element with the current material number Care should be taken to minimize the number of history variables since for very large problems the memory requirements can bec
29. ative along the axis of the element and SHP 2 i the shape function N In calculations integrals are represented as 1 EN Nids f ANO NaO XJACE de 6 5 L 1 and quadrature may be used for evaluation Calculation of natural coordinate derivatives only may be obtained with the call CALL SHAP1DN S SHP NEL where Parameter Description 5 Natural coordinate SHP 2 NEL Shape function and derivative NEL Number element nodes 2 or 3 where SHP 1 i contains N g and SHP 2 i the shape function N Cubic Hermitian interpolation e g for use in straight linear beam elements given by w NY n NY z NEO N 0 6 6 is obtained using the call CHAPTER 6 UTILITY ROUTINES 68 CALL SHP1DH S LEN SHPW SHPT where Parameter Description S Natural coordinate LEN Length of the element 2 node SHPW 4 2 Shape functions for z SHPT 4 2 Shape functions for 6 The arrays are evaluated as follows 1 SHPW 1 i SHPT 1 i are first derivatives e g Ni 2 SHPW 2 i SHPT 2 i are second derivatives e g Ni xs 3 SHPW 3 i SHPT 3 i are third derivatives e g Nj srx and 4 SHPW 4 i SHPT 4 i are shape functions e g N 6 2 2 Shape functions in two dimensions Two dimensional Co isoparametric interpolation on quadrilaterals of linear quadratic and cubic order may be obtained using the subprogram call CALL SHP2D SS XL SHP XJAC NDM NEL IX
30. ay its values are set to zero Thus if the array is to be used only once it need not be set to zero before accumulating additional values If the array is to be reused or resized see below it must be reinitialized prior to accumulating any additional values Use of these subprograms controls the assignment of memory space for all arrays such that no CHAPTER 3 ALLOCATING ARRAYS 10 conflicts occur between hr and mr referenced arrays Each routine which makes direct reference to an allocated array using a pointer e g hr np 43 or mr up 1 must contain include files as include pointer h include comblk h As an example for the use of the above allocation scheme consider a case where it is desired to allocate a real double precision array with length NUMNP number of nodes in mesh and an integer array with length NUMEL number of elements in mesh The parameters NUMNP and NUMEL are contained in COMMON CDATA and available using the include file cdata h The new arrays aredefined using the temporary names TEMP1 and TEMP2 which have numerical locations 111 and 112 respectively The two arrays are allocated using the statements setvar setvar palloc 111 1 numnp 2 palloc 112 TEMP2 numel 1 where the last entry indicates whether the array is REAL 8 2 or INTEGER 1 These arrays are now available in any subprogram by specifying the pointer h and comblk h include files and referencing
31. ble UCT logical uprt real 8 ct1 3 C Set command word if bcomp uct p1t0 4 then uct else c User plot command statements are placed here endif end The parameters CTL 3 are used to pass the three parameter values read respectively Again the name xxxx should be selected to not conflict with existing plot command names and will appear whenever HELP is issued Two plot utilities are available for placing lines on the screen These are named DPLOT and PLOTL The calling form for DPLOT is given as call dplot s1 s2 ipen where s1 s2 are screen coordinates ranging from 0 to 1 Similarly the calling sequence for PLOTL is CHAPTER 4 USER FUNCTIONS 26 Number Color 0 Black 1 White 2 Red 3 Green 4 Blue 5 Yellow 6 Cyan 7 Magenta 8 Orange 9 Coral 10 Green Yellow 11 Wheat 12 Royal Blue 13 Purple 14 Aquamarine 15 Violet Red 16 Dark Slate Blue 17 Gray 18 Light Gray Table 4 1 Color Table for Plots call plotl x1 x2 x3 ipen where x1 x2 x3 are coordinates values of the mesh The value of ipen ranges from 1 to 3 1 starts a filled panel 2 draws a line from the current previous point to the new point 3 moves to the new point without drawing a line If a filled panel is started it must be closed by inserting the statement call clpan Lines are drawn or panels filled in the current color A color is set using the statement call
32. contour values The known types of plots are 1 Point element with one node obtained by call CALL PLTPT1 IEL 2 Line element with two nodes obtained by call CALL PLTLN2 IEL and for three node elements CALL PLTLN3 IEL 3 Triangular element with 3 nodes obtained by call CALL PLTRI3 IEL and for 6 nodes obtained by call CALL PLTRI6 IEL CHAPTER 6 UTILITY ROUTINES 73 4 Quadrilateral element with 4 nodes obtained by call CALL PLQUD4 TEL for 8 or 9 node elements the plot call is CALL PLQUD8 TEL and for 12 or 16 node quadrilaterals the call is CALL PLTQ16 IEL 5 Tetrahedral element with 4 nodes obtained by call CALL PLTET4 IEL and for 10 node tetrahedra the call is CALL PLTET10 TEL 6 Brick element with 8 nodes obtained by call CALL PLBRK8 IEL and for 20 or 27 node bricks the call is CALL PLBRK27 IEL Using these and internal extraction of element surfaces the program is able to make some hidden surface plots in three dimensions 6 4 2 Element data plots Users may construct plots within their elements i e an ELMTnn and access using the plot command PLOT PELE v1 v2 v3 In interactive mode in the plot environment it is only necessary to enter CHAPTER 6 UTILITY ROUTINES 74 PELE vi v2 v3 The values entered in v1 v2 v3 are optional and are passed to the element through a common block as REAL 8 ELPLT COMMON ELPDAT ELPLT 3 The PELE option ca
33. d in s for the various cases a parameter array ctan 3 is passed to the subprogram in labeled common eltran When the task parameter isw is 3 the values in the ctan array are interpreted according to Table 5 9 Thus in solid mechanics applications the tangent matrix is defined in an element routine as S ctan 1 K ctan 2 C ctan 3 M where K is the stiffness matrix C is the damping matrix and M is the mass matrix For non linear applications these matrices normally are computed with respect to the current values of the available solution parameters The values provided in the ctan array are set by FEAP according to the active transient solution option For a static option both ctan 2 and ctan 3 are zero For options integrating first order differen tial equations in time only ctan 3 will be zero For options integrating second order differential equations in time all the parameters are non zero 5 3 Example 2 Node Truss Element An element routine carries out tasks according to the value assigned to the parameter isw as indicated in Table 5 2 To describe basic steps to program the various tasks CHAPTER 5 ADDING ELEMENTS 42 defined by isw we consider next the problem of a 2 node linear elastic truss element for small deformation applications The element is described in sufficient generality to permit solution of both two and three dimensional truss problems 5 3 1 Theory for a Truss The governing equations for a typical truss
34. d in the subprogram PALLOC The subroutine PGETD also may be used to retrieve internal data arrays by NAME for use in user developed modules For example if a development requires the nodal coordinate data the call integer xpoint xlen xpre logical flag call pgetd X xpoint xlen xpre flag will return the first word address in memory for the coordinates as xpoint the length of the array as xlen and the precision of the array as xpre If the retrieval is successful flag is returned as true whereas if the array is not found it is false The precision will be either one 1 or two 2 for INTEGER or double precision REAL 8 quantities respectively Thus the above coordinate call will return xpre as 2 and xlen will be the product of the space dimension of the mesh and the total number of nodes in the mesh The first coordinate z may be given as x1 hr xpoint any other coordinates at nodes may also be recovered by a correct positioning in later words of hr For example y is located at hr xpoint 1 The use of pgetd can lead CHAPTER 3 ALLOCATING ARRAYS 14 to errors for situations in which the length of arrays changes during execution since in these cases the value of the pointer xpoint can change For such cases a call to pgetd must be made prior to each reference involving xpoint On the other hand reference using the pointers defined in arrays NP or UP are adjusted each time an array changes size However users must
35. divided into two 4 entry quantities if required Indeed the space may be used in anyway necessary provided no conflict in the way each time value is associated to the data Note that sufficient storage is available to define integration methods for which the rotation is defined at an intermediate time tnya For a typical node n in the mesh the location of the entries in the x1g array are obtained from ng mropt n 2 and the routine UROTmm is called as call urotmm xlg 1 1 ng xlg 1 2 ng xlg 1 3 ng amp xlg 1 4 ng xlg 1 5 ng amp xlg 4 4 ng xlg 4 5 ng amp xlg 7 4 ng xlg 7 5 ng du j n tsu where du j n are the solution increments from the solver and tsw is the time update switch which is set according to tsw 1 Initialize for new time step tsw 2 Update within a time step tsw 3 Back up to beginning of time step The entry u j n is the location for the first entry in the rotation vector 6 CHAPTER 5 ADDING ELEMENTS 53 5 7 2 Local nodal rotation treatment When each element that is associated with nodal rotation parameters is called the rotation data is transferred to local storage in a manner similar to treatment of trans lations The local data is passed to each element using the common erotas defined above The entries in the local arrays are extracted from the global array according to xln 1 9 n1 1 xlg 1 9 1 ng xln 1 9 n1 2 xlg 1 9 2 ng x1n 1 9 n1 3 xlg 1 9 3 ng xln 1 9 n1 4 x
36. e nodal connection array but is used to return the list of unused degrees of freedom Utility routines are also provided to assist users in providing the necessary list of nodes needed to properly draw the mesh each element type during plot outputs The names of the routines are listed in Table 5 10 and each routine is called as call plname iel where iel is an integer parameter defined in common eldata If no call to a subpro gram is included each element is assumed to be a 4 to 9 node quadrilateral and default drawing order will be assigned Users may construct their own drawing order also by following the steps employed in one of the routines defined in Table 5 10 5 4 2 Task 6 Options The TPLOt solution command includes an option to save specific element quantities e g stress strain etc This option is implemented for user elements by including the common CHAPTER 5 ADDING ELEMENTS 47 Routine Name Description PLTLN2 2 node line element PLTRI3 3 node triangular element PLQUD4 4 node quadrilateral element PLTRI6 6 node triangular element PLTET4 4 node tetrahedron element PLBRK8 8 node brick element Table 5 10 Element Plot Definition Subprograms real 8 tt common elplot tt 100 and then inserting the statement tt i value at an appropriate location in the isw 3 task For example if it were desired to save the force and strain in the truss element the statements tt 1 EAteps Elem
37. ee we e Kee ae ah 22 5 1 FEAP Element Subprogram Parli aa arta 27 5 1 FEAP Element Subprogram 28 5 2 FEAP Element Common Blocks waa pag tere teren dte Here BO 57 5 3 FEAP Element Common Blocks using Includes 58 5 4 2 Node Truss Element en Saer de S M x B 58 5 5 Element residual for two node truss 59 5 6 Truss Tangent Matrix Part 1 anar an ne BA ee 60 5 7 Truss Tangent Matrix 0330000000 61 5 8 Element variable projection routine 62 iii List of Tables 3 1 3 2 3 3 3 4 3 5 4 1 5 1 5 2 5 3 5 4 5 9 5 6 5 7 5 8 5 9 5 10 5 11 6 1 6 2 6 3 6 4 Mesh Array Names Numbers and Sizes 11 Solution Array Names Numbers and Sized 11 Element Array Names Numbers and Sizes 11 Element connection array IX use for element e 12 Element control array IE use for material number ma 12 Color Tablefor Plots a dos a GA ae et ae Se Ma aa AN ag 26 Arguments of FEAP Element Subprogram 29 Task Options for FEAP Element Subprogram R Required O Optional E For eigensolutions 30 FEAP common block definitions Waa 32 FEAP common block definitions 33 Material Parameters ip Leida Pa eb Sie Ra ee 34 Material Parameters sate bee s ned aloni 35 Material Parameters era beta
38. ement a Newton method it is necessary to linearize the residual equation For FEAP the Newton equation may be written as OR RT R da 0 i i Oa j where a is one of the variables at time t 1 e g Us tn 1 We define s R w k da and solve sda RP The solution is updated using k 1 k k a da In the above k is the iteration number for the Newton algorithm To start the solution for each step FEAP sets at alta where a quantity without the k superscript represents a converged value For a linear problem Newton s method converges in one iteration Computing the residual after one iteration must yield a zero value to within the roundoff of the computer used For non linear problems a properly implemented Newton s method must exhibit a quadratic asymptotic rate of convergence Failure of the above performance for linear and non linear cases implies a programming error in an implementation or lack of a consistently linearized algorithm i e S is not an exact derivative of the residual In a non linear problem Newmark s method may be parameterized in terms of in crements of displacement velocity or acceleration From the Newmark formulas the relations du BAT da and dv yAtda CHAPTER 5 ADDING ELEMENTS 39 define the relationships between the increments Note that only scalar multipliers involving 0 y and At are involved between the different measures The ta
39. ent axial force tt 2 eps Element axial strain could be placed anywhere after the stress and strain are defined These values would be output by using a solution command sequence such as batch tplot end stress nn i saves force for element nn stress nn 2 saves strain for element nn show writes tplot items to output file 5 5 Projection of element variables to nodes The STREss NODE solution command and the PLOT STRE command require a projection of element variables to nodes A continuous stress field is assumed to obtain the nodal CHAPTER 5 ADDING ELEMENTS 48 values Accordingly o Na Ta where is any value which is to be projected to nodes e g a stress or strain Na are shape functions for the element type considered and ca nodal values of the projected quantity The projection routine uses a diagonal weight matrix to project the values For simple elements the matrix is computed by a procedure identical to mass lumping For example Moa i Na dQ Q defines a row sum form of projection 3 Using the above results in the set of equa tions and a least square fit with the finite element values 6 gives the equation set Ma Ga aj Nga dn Q This defines nodal values for projected quantities Since the coefficient matrix is di agonal the solution to the set of equations for each component is trivial The actual solution is performed automatically by FEAP To permit each element to project its own qu
40. eration in direction x v e b is a loading in direction x per unit length and e de is a virtual strain along the truss axis For a straight truss member the displacement along the axis us may be expressed in terms of the components in the directions x as d us l u s t Y liui s t i 1 where t is time u is the displacement vector with components u l is a unit vector along the axis of the member with direction cosines l defined by Ox _ vig Ti p s L d Drz Y zo Ta i 1 and zi tj are the coordinates of nodes 1 and 2 respectively The displacement components are interpolated on the 2 node truss member as uls 1 ualt uo t 6 7 CHAPTER 5 ADDING ELEMENTS 44 in which uj ujo are the displacements at nodes 1 and 2 The virtual displacements are obtained from the above by replacing u by du etc The truss strain is I x Os Using the interpolations for the displacement components yields 1 d Ess T2 Ax Au i 1 where Az Li Lil L L and Au Uig Ujj Thus in matrix form the strain is 1 i 2 Ess d i2 Az Ax 7 1 Using the above displacement interpolations the variational equation for the truss may be expressed in matrix form as BIL un duo I 7 a oss Ads ee pA 1 ds pa AS 2 R balam FEAP constructs the finite element arrays from the element residual
41. gration compute shape functions and their derivatives etc 6 1 Numerical quadrature routines Details on quadrature formula types and the layout and location of points and weights may be found in standard references 5 6 2 3 Here only the description of subroutine calls is included together with the available options on number of points 6 1 1 One dimensional quadrature Line integrals may be evaluated using Gaussian quadrature in which the approximation to an integral is given as 1 L FO de Y 6 W 6 1 1 1 1 where amp are quadrature points and W are the weights to be applied at each point The weights satisfy the condition S Wi 2 6 2 CHAPTER 6 UTILITY ROUTINES 64 The Gauss Legendre formula has points amp which are all less than unity The subpro gram call CALL INTID L SG WG in which L is assigned an integer value between 1 and 5 returns the points in the array SG and the weights in WG both of which are of type REAL 8 The Gauss Legendre formula integrates exactly polynomials up to order 2 L 1 The Gauss Lobato formula has two of its points at 1 and 1 with the remainder in the interior of the interval A routine to perform quadrature is obtain by using the call CALL INTIDL L SW in which L is assigned an integer value between 1 and 6 The values of the points and weights are returned in the two dimensional array SW Points in SW 1 and weights in SW 2 6 1 2 Two di
42. hape type 1 rectangles 2 tube 3 Wide flange 4 Channel 5 Angle 5 Circle 101 126 Shape data 127 Surface convection h 128 Free stream temperature Tx 129 Reference absolute temperature 130 nseg Number of hardening segments 131 148 Segment data sets epYjsoHkin 149 Total variables on frame section 150 Piezoelectric flag 151 159 Piezoelectric data 160 Initial stress flag 161 166 0 Initial stresses constant 167 Tension compression only indicator 170 C Fung pseudo elastic model modulus 171 41 Fung model energy parameter 172 ag Fung model energy parameter 173 a3 Fung model energy parameter 174 ag Fung model energy parameter 175 as Fung model energy parameter 176 ag Fung model energy parameter 177 az Fung model energy parameter 178 ag Fung model energy parameter 179 ag Fung model energy parameter 180 181 Viscoplastic rate parameters 182 Nodal quadrature parameters 183 Bm Mz Mc mass scaling factor 184 ic Estimate on maximum wave speed 185 Augmentation switch 187 Implicit 0 Explicit 1 element integration 188 200 Unused Table 5 8 Material Parameters CHAPTER 5 ADDING ELEMENTS 38 and Vi tn4 1 Vilta At 1 y ai tn vai tn41 with u v and a being the set of displacements velocities and accelerations at node 1 respectively A Newton method is commonly adopted to solve a non linear or linear problem To impl
43. hapter 3 can be used to develop user routines for manipulation of the mesh data For example if a user has added the specification of information by coordinates it may later be necessary to associate the data with specific node numbers This can be accomplished using a manipulation command which searches for the node number whose coordinates are closest to the specified location 4 4 Solution Command Functions UMACRn In a similar manner users may add solution commands to the program by adding a routine with the name UMACRn where n ranges from 0 to 9 subroutine umacrO lct ctl uprt User solution command function implicit none include umaci h I Contains the variable UCT logical uprt character lct 15 real 8 ct1 3 C Set command word if pcomp uct mac0 4 then uct else c User command statements are placed here endif end CHAPTER 4 USER FUNCTIONS 25 The parameters LCT and CTL are used to pass the second word of a solution command and the three parameter values read respectively Again the name xxxx should be selected to not conflict with existing solution command names and will appear whenever HELP is issued 4 5 Plot Command Functions UPLOTn In a similar manner users may add new plot commands to the program by adding a routine with the name UPLOTn where n ranges from 0 to 9 subroutine uplotO ctl uprt c User plot command function implicit none include umac1 h Contains the varia
44. he eleventh parameter of the ELMTnn subprogram A short description of the task carried out by each value as currently implemented is shown in Table 5 2 To use the basic features available in FEAP it is necessary to program tasks labeled as R shown above If elements have local variables that need to be retained between subsequent time steps history variables may be defined as described in Section 5 6 In this case it is necessary to code task 12 if special transformations of the variables are required otherwise merely return with no changes and if any of the parameters have non zero initial values task 14 is used to set these values zero values are set by default Finally if special plotting options are desired it may be necessary to program task 20 note that contours for element variables such as stress strain etc are developed from task 8 Parameter Description d Element data parameters Moduli body loads etc ul ndf nen j Element nodal solution parameters nen is number of nodes on an element max j 1 Displacement ub j 2 Increment ul Un j 3 Increment u ur j 4 Rate oP j 5 Rate a 6 Rate vn xl ndm nen Element nodal reference coordinates ix nen Element global node numbers tl nen Element nodal temperature values s nst nst Element matrix e g stiffness mass r ndf nen Element vector e g residual mass may also be used as r nst ndf Number unknowns max pe
45. heck element for errors Negative jacobian etc a elseif isw eq 3 then Return Element coefficient matrix and residual s nst nst element coefficient matrix r ndf nen element residual hr nh1 history data base previous time step hr nh2 history data base current time step hr nh3 history data base time independent 000000 elseif isw eq 4 then c Output element quantities e g stresses elseif isw eq 5 then c Return Element mass matrix c s nst nst consistent matrix r ndf nen diagonal matrix elseif isw eq 6 then c Compute residual only c r ndf nen element residual elseif isw eq 7 then c Return Surface loading for element c s nst nst coefficient matrix c r ndf nst nodal forces elseif isw eq 8 then c Compute stress projections to nodes diagonal p nen projection weight wt nen c s nen projection values st nen default project 8 quantities endif end Figure 5 1 FEAP Element Subprogram Part 2 11 items which are briefly described in Table 5 11 Note in Table 5 1 that FEAP transfers the values for most of the solution parameters in array UL NDF NEN at time tn a Where a denotes a value between 0 and 1 The value of a is 1 ie values are reported for time t 1 unless generalized midpoint integration methods are used For the present we will assume a is 1 CHAPTER 5 ADDING ELEMENTS 29 FEAP carries out tasks according to the value of the parameter ISW passed as t
46. hould be double preci sion and IPR is set to 2 For systems in which REAL 8 variables are single precision with the same working length as integer variables the IPR parameter is set to 1 Any error in setting this parameter may lead to incorrect behavior of the program consequently do not reset the parameter to single precision unless a careful assessment of compiler behavior for REAL 8 variables has been made By placing an alphanumeric version of each manual page in a separate file which has the name of the command and a t extender e g coor t for the mesh coordinate input command it is possible to read each page during execution using the HELP name command where name is the command name whose manual page is to be read For CHAPTER 1 INTRODUCTION 2 this option to work properly it is necessary to define the path name to each manual page in the FEAP82 module For example file 1 c Feap Manual Mesh file 2 c Feap Manual Macr file 3 c Feap Manual Plot defines a typical path for a PC system Each system requires a proper path definition FEAP will add the requested command name to each of the above paths to find mesh solution or plot commands Normally FEAP reads each input data line as text data and checks each character for the presence of parameters expressions and constants For very large data sets this parsing of each instruction can consume several seconds of compute time If all data is norma
47. i 25 doz Material moduli 26 d31 Material moduli 27 912 Material moduli 28 923 Material moduli 29 931 Material moduli 30 htype Heat flag Table 5 5 Material Parameters 34 CHAPTER 5 ADDING ELEMENTS Parameter Name Description 31 4 Orthotropic angle x principal axis 1 32 A Area cross section beam truss 33 T Inertia cross section beam truss 34 Iz Inertia cross section beam truss 35 h Inertia cross section beam truss 80 J Polar inertia cross section beam truss T Ki Shear factor plate 38 Ko Shear factor plate 39 Non linear flag beam truss 40 Inelastic material model type 41 Yo Initial yield stress Mises 42 Wns Final yield stress Mises 43 G Exponential hardening rate 44 Hiso Isotropic hardening modulus linear 45 Hun Kinematic hardening modulus linear 46 Yield flag 47 Bi Orthotropic thermal stress 48 z Orthotropic thermal stress 49 B3 Orthotropic thermal stress 50 Error estimator parameter l yy Viscoelastic shear parameter 52 Ti Viscoelastic relaxation time 53 va Viscoelastic shear parameter 54 T Viscoelastic relaxation time 55 13 Viscoelastic shear parameter 56 73 Viscoelastic relaxation time 57 nvis Number of viscoelastic terms 1 3 58 Damage limit 59 Damage rate 60 k Penalty parameter Table 5 6 Material Parameters 39 CHAPTER 5 ADDING ELEMENTS Parameter
48. ites to output file c and or write LABEL list writes to screen endif debug This device supplements use of available debuggers on the computer 1 2 Uses of Common and Include Statements FEAP contains many COMMON statements that are used to pass parameters and small array values between subprograms For example access to the debugging parameter debug is facilitated through common debugs Users may either place the common statement as well as data typing statements directly in the routine or may use an include statement For debugging the statement would be CHAPTER 1 INTRODUCTION 4 include debugs h which during compilation would direct the precompiler to load the current common statement from this file In FEAP all include files have the same name as the common with an added extender h The only exception is a dummy blank common which uses the file name comblk h which is defined as real 8 hr integer mr common hr 1 mr 1 The dummy arrays hr 1 and mr 1 serve to pass all dynamically allocated arrays be tween subprograms using a pointer array contained in the common array named or for user defined arrays in up located in the include file pointer h See Section 3 for more details on use of pointers All include files are located in the directories include and include integer4 or include integer8 The subdirectory integer4 defines 32 bit integer parameters and integer8 defines 64 bit integer
49. lg 1 9 6 ng rots 1 3 n1 1 xlg 1 3 4 ng rots 1 3 n1 2 x1g 1 3 5 ng rvel 1 3 n1 1 xlg 4 6 4 ng rvel 1 3 n1 2 xlg 4 6 5 ng racc 1 3 n1 1 xlg 7 9 4 ng racc 1 3 n1 2 xlg 7 9 5 ng where nl is a local node number between 1 and 9 the maximum provided in the current erotas and ng is the global node number associated with each local number Using the above data structure one can program the updates in any manner that does not conflict with the time treatment The only interface to FEAP is through how the increment du 4 5 n is defined 5 8 Energy Computation FEAP elements provide an option to accumulate the total momenta and energy during the solution process The values are accumulated in the array EPL 20 when the switch parameter isw is 13 and written to a file named Pxxxx ene where xxxx is extracted from the problem input filename whenever the solution command TIME is used The array EPL 2 is in the common block named ptdat6 which has the structure 1 8 1 integer iepl neplts common ptdat6 ep1 20 iep1 2 20 neplts For problems in solid mechanics the linear momenta are stored as follows CHAPTER 5 ADDING ELEMENTS 54 Component Description EPL 1 EPL 3 Linear momenta EPL 4 EPL 6 Angular momenta EPL 7 Kinetic energy EPL 8 Stored energy EPL 9 Work by external loads EPL 10 Total energy Table 5 11 Momenta and Energy Assignments The linear momenta are comp
50. lls each element with the switch parameter ISW 20 Users merely code whatever option they wish to include within their element module The standard color table is available through use of the subroutine call CALL PPPCOL ICOL 0 in which ICOL designates the color to be assigned according to Table 6 3 An exception occurs for PostScript outputs where black and white are switched since the background then is assumed to be white ICOL COLOR ICOL COLOR 0 Black 10 Green yellow 1 White 11 Wheat 2 Red 12 Royal blue 3 Green 13 Purple 4 Blue 14 Aquamarine 5 Yellow 15 Violet red 6 Cyan 16 Dark slate blue 7 Magenta 17 Grey 8 Orange 18 Light grey 9 Coral Table 6 3 Color pallet for FEAP plots A straight line segment may be drawn to the screen in the current color between the coordinates 1 yn zi and z Ya 22 using the commands CALL PLOTL X1 Y1 Z1 3 CALL PLOTL X1 Y1 Z1 2 Here the basic command is CHAPTER 6 UTILITY ROUTINES 75 IP Action 1 Start panel fill 2 Move to point 3 Draw to point Table 6 4 Values for control of plots CALL PLOTL X1 Y1 Z1 IP where the three cartesian coordinates relate to mesh coordinates not screen values and IP is a parameter defined according to Table 6 4 The perimeter of a panel is drawn with standard line drawing commands starting with CALE PLOTL X1 Y1 Z1 1 and continuing with a seque
51. lly provided as numerical data without use of any parameters or expressions the input time may be reduced by setting the value of the logical variable COFLG in FEAP82 to false FEAP will automatically switch to parsing mode if any record contains non numerical data item It is also possible to use the PARSe and NOPArse commands to set the appropriate mode of data input During the input of plot commands FEAP has the option to either set input options automatically DEFAult mode or to read the values or range of contours to plot The default mode of operation may be assigned in the FEAP82 module by setting the variables DEFALT and PROMPT Setting DEFALT to true indicates that all default options are to be set automatically If DEFALT is set false a prompt for contour intervals may be requested by setting PROMPT to true FEAP has options to produce encapsulated PostScript output files in either gray scale of in color The default mode may be established by setting the variable PSCOLR and PSREVS Setting PSCOLR true indicates the PostScript files will be in color unless set otherwise by the PLOT COLOr data command The PSREVS variable reverses the color sequence Arrays in FEAP are dynamically allocated during execution Thus it is possible to define and destroy arrays as well as to increase or decrease the size of an array A pa rameter is provided to control when an array is to be decreased in size The parameter is INCRED and an array is decrea
52. m INPT2D This subroutine is accessed by the statement CALL INPT2D D TDOF NEV TYPE where D is the array storing the material parameters TDOF is returned as the parameter to access temperature NEV is the number of element history variables to allocate to NH1 and TYPE is the element type This routine inputs the commands as described in the user manual and stores the data for each material set into the D array elements as described in the following tables CHAPTER 5 ADDING ELEMENTS 32 Variable Definition o Page eject option head Title record numnp Number of mesh nodes numel Number of mesh elements nummat Number of material sets nen Maximum nodes element neq Number active equations ipr Real variable precision nstep Total number of time steps niter Number of iterations current step naugm Number of augments current step titer Total iterations taubm Total augments iaugm Augmenting counter iform Number residuals in line search dm Element proportional load n Current element number ma Current element material set mct Print counter iel User element number nel Number nodes on current element tt Element stress values for TPLOt bpr Principal stretch ctan Element multipliers ut Element user values for TPLOt Table 5 3 FEAP common block definitions 5 2 Non linear Transient Solution Forms Before describing the steps in developing an element we summarize first the basic struct
53. mensional quadrature Two dimensional quadrature on quadrilateral domains may be performed by repeated one dimensional integration The two dimensional integrations are approximated by remaean fem W 6 3 where L is the total of all quadrature points A routine to compute n x n order Gauss Legendre quadrature is obtained by the call CALL INT2D L LINT SW where L is assigned to the number of points in each direction LINT is returned as the total number of points and SW 3 is an array containing the points and weights according to SW 1 1 contains values of the points SW 2 1 contains values of the points 7 and SW 3 1 contains values of the weights W Two dimensional quadrature on triangles may be performed using the subprograms call CHAPTER 6 UTILITY ROUTINES 65 Type Number Location Points 1 1 Centroid O h 3 3 Mid sides O h 3 3 Interior O h 4 4 Interior O h Negative Wt 6 6 Nodal O h 6 6 Interior O h rd 7 Interior O h5 7 7 Noda O h9 12 12 Interior O h 13 13 Interior O h3 Negative Wt Table 6 1 Quadrature for triangles CALL TINT2D L LINT SW where L is a type indicator LINT returns the number of points and SW 4 is an array which returns three area coordinates and the quadrature weight SW 1 1 returns the area coordinate Li as defined in 2 31 SW 2 1 returns the area coordinate L a SW 3 1 retur
54. n solution Element residual array storage 40 Element tangent array storage 39 Non linear problems 38 Tangent array definition 39 Numerical integration 63 One dimension 63 Three dimension 65 Two dimension 64 Plot command UPLOT subprogram 76 UPLOTn subprogram 25 Plot utility function Fills 75 Lines 74 Plot utility functions Color 26 Fills 25 Lines 25 Precision Integer 1 Real 1 Principal stress Three dimensional 71 Two dimensional 71 Problem size 1 Quadrature 63 One dimension 63 Three dimension 65 Two dimension 64 Setting colors 26 Setting options Help level 2 Parsing statements 2 Plot prompts 2 PostScript 2 Shape functions 66 79 Brick isoparametric 69 One dimensional cubic Hermitian 67 One dimensional isoparametric 67 One dimensional natural derivative 67 Quadrilateral isoparametric 68 Tetrahedral isoparametric 70 Three dimensional isoparametric 69 Trianglular isoparametric 69 Two dimensional isoparametric 68 Solution array names 11 Solution command UMACRn subprogram 24 Table of colors 26 Transient solution CTAN definition 41 Newton method 38 Non linear 32 Truss element Newmark method 45 Non linear theory 54 Tangent modulus 56 Theory 41 Variational equation 44 Weak form 44
55. nce of draw commands CALL PLOTL Xi Yi Zi 2 however no lines appear on the screen and the fill of each panel is completed by the statement CALL CLPAN It should be noted that all plots within FEAP are performed in three dimensions For two dimensional problems no z coordinates are available in the XL NDM NEN array and hence it is necessary to assign zero values for the z coordinates before calling a plot subprogram If a perspective view has been requested a full use of a i Yi zi Specification is made In this case a user may wish to pass the value of some solution variable as the z value scaled so that it will make sense relative to the x y coordinate values Similarly if deformed plots are being performed it is necessary to add scaled displacements to the coordinates The current value of the scaling parameter i e variable CS is available in labeled common PVIEW In this case one can add the statements assuming here that the displacements correspond to the coordinate directions DO NE 1 NEL CHAPTER 6 UTILITY ROUTINES 76 DO I 1 NDM XP I NE XL I NE CS UL I NE END DO I END DO NE NEL is the number of connected nodes to each element and is passed through labeled common ELDATA before performing any deformed plots and then plot the appropriate values of XP Indeed this may always be performed as the value of CS will be zero for an undeformed plot 6 4 3 Other user plots It is also possible fo
56. ndf the marimum number of degree of freedoms at any node in the problem and nen the maximum number of nodes on any element The ordering of the unknowns into nst must be carefully aligned in order for FEAP to properly assemble each element matrix into the global tangent The ordering is such that sub matrices are defined for each node attached to the element Thus Su Si Sis Sa 922 S23 Sa S32 Sas CHAPTER 5 ADDING ELEMENTS 40 where S is the sub matrix for nodal pairs 7 7 Each of the sub matrices is a square matrix of the size of the maximum number of degree of freedoms in the problem which is passed to the subprogram as ndf Thus Si Si Sis S Sa Sap S33 Si S32 933 Sp nap in which Si is an array coefficient for nodal pair i j for the degree of freedom pair a b In FEAP the element residual may be stored as one dimensional array which is dimen sioned r nst with entries stored in the same order as the rows of the element tangent matrix or as a two dimensional array which is dimensioned as r ndf nen The one dimensional form of the residual is given as Ri R R R where the entries in each submatrix are given as Ri FR R fg Ria The two dimensional form places the entries R as columns Accordingly R R R Ra The two forms for defining the residual r are equivalent based on the Fortran ordering of information into double subscript arrays If ndf is larger than needed for the element
57. ngent matrix for the transient problem using Newmark s method may be ex pressed in terms of the incremental displacement velocity or acceleration As an example consider the case where the solution is parameterized in terms of increments of the displacements i e a is the displacement vector For this case the tangent matrix is we do not show dependence on the iteration k for simplicity of notation Sdu Piga OR 976 OR Oa 7 Ov du Ou 7 Note that from the Newmark formulas dar 1 5 e _ eda _ Y 5 du GARE du da du in which is the Kronnecker delta identity matrix for the k j nodal pair From the residual we observe that Cy My du 2 Ov i da define the tangent stiffness damping and mass respectively Thus for the Newmark algorithm the total tangent matrix in terms of the incremental displacements is Y 1 B At 4 For other choices of increments the tangent may be written in the general form Sij Ki 4 Sis Ci K C2 C3 Mi where the c are scalar quantities involving the integration parameters of the method selected and At Thus any one step integrator may be considered and will affect only the specification of the constants in the tangent m FEAP the element tangent matrix S is stored as a two dimensional array which is dimensioned as s nst nst where nst is the product of ndf and nen with
58. ns the area coordinate Lg SW 4 1 returns the weight Wy Table 6 1 describes the admissible types number and location of quadrature points 6 1 3 Three dimensional quadrature Three dimensional quadrature on brick domains may be performed by repeated one dimensional integration The three dimensional integrations are approximated by 1 L ein O dI YO FE 6 4 m l 1 where L is the total of all quadrature points A routine to compute n x n x n order Gauss Legendre quadrature is obtained by the call CALL INT3D L LINT SW where L is assigned to the number of points in each direction LINT is returned as the total number of points and SW 4 is an array containing the points and weights according to SW 1 1 contains values of the points SW 2 1 contains values of the CHAPTER 6 UTILITY ROUTINES 66 Type Number Location Points 1 1 Centroid O h 1 4 Nodal h 2 4 Interior O h 3 5 Interior O h Negative wt 4 11 Interior O h Negative wt 4 11 Nodal O h3 5 14 Interior O h 6 16 Interior O h 8 8 Node amp Face O h Table 6 2 Quadrature for tetrahedra points 7 and SW 3 1 contains values of the points and SW 4 1 contains values of the weights W Three dimensional quadrature on tetrahedra may be performed using the subprograms call CALL TINT3D L LINT SW where L is a type indicator LINT returns the number of p
59. ntities the virtual strain may be written as ddu a Na g 1 ss where After differentiation of the displacement field the discrete form of the virtual strain is given by 1 7 du du E CHAPTER 5 ADDING ELEMENTS 56 Substituting the above virtual strain expression into the weak form gives the modified residual expression _1 1 Yi L 2 1 s x n m IE The tangent tensor is obtained by linearizing the residual as shown previously The only part which is different is the term with o Noting that Odu d ss DE and dou n n n n ddu Os i Os If the n are constructed as column vectors then the tensor product becomes a matrix defined as d ess T T G n Qn Ono nini nons With these definitions the tangent matrix for the non linear problem is given as EA Gi Oss A Gi ad E Aa 9 gi 2 J Gij Gi Notice that for the linear problem thus the only difference between the linear and non linear problem is the definition of ss in terms of displacements the modification for geometric effects for the g and the second term on the tangent matrix which is sometimes called the geometric stiffness part CHAPTER 5 ADDING ELEMENTS 57 character 4 common bdata integer common cdata integer common counts integer common counts real 8 integer common eldata real 8 common elplot
60. ny of the hr nh1 hr nh2 or hr nh3 types the new quantities will be returned to the H history for the element only for isw tasks where residuals are being formed for a solution step i e solution command FORM TANG 1 or UTAN 1 and for history reinitialization during a time update i e solution command TIME These access the task options isw equal to 3 or 6 and 14 respectively If a user adds a new option for which it is desired to save the history variables it is necessary to set the variables hflgu and h3flgu to true as required if no update is wanted the variables should be set to false These parameters are located in logical hflgu h3flgu common hdatam hflgu h3flgu 5 7 Elements with Finite Rotation Parameters When considering structural elements that undergo large displacements it is usually necessary to treat the rotation parameters for large angle changes The nodal pa rameters for this case are a rotation vector 0 and the finite rotations are given as an CHAPTER 5 ADDING ELEMENTS 51 orthogonal matrix A A Anti exp 6 An in which 6 denotes a skew matrix given as 9 0 6 0 0 A 0 A 0 The actual method used to update the rotations and their increments must be specified when writing the element module ELMTnn and is performed by a user subprogram named UROTmm where mm is a number between 01 and 10 To specify which routine is to be used it is necessary to include the statement rotyp mm
61. o input and the nn specifies the number of real data values to input The value for parameter nt or nn may be zero Thus the use of errck tinput text 0 td nn CHAPTER 2 DATA INPUT AND OUTPUT 7 is equivalent to errck pinput td nn Text variables may be converted to numerical REAL 8 form using the subroutine call call setval text nc td where text is a string with nc characters and td a REAL 8 variable The text string can contain any parameters expressions or numerical constants which evaluate to a single value 2 2 Array Outputs Two subprograms exist to output arrays of integer and real double precision data The routine MPRINT is used to output real data and is accessed by the statement call mprint array nrow ncol ndim label where array is the name of the array to print nrow and ncol are the number of rows and columns to output ndim is the first dimension on the array and label is a character label which is added to the output For example the statements realx8 aa 8 6 call mprint aa 2 4 2 3 8 AA outputs a 2 x 3 submatrix from the array aa starting with the entry aa 2 4 The output entries will be ordered as the terms aa 2 4 aa 2 5 aa 2 6 aa 3 4 aa 3 5 aa 3 6 The MPRINT routine adds row and column labels as well as the character label The routine IPRINT is used to output integer data and is accessed by the statement call iprint array nrow ncol ndim label where all p
62. obian transformation from xyz to n SHP 4 4 Shape functions and derivatives The array SHP stores the values in the same order as for the brick element 6 3 Eigenvalues for 3 x 3 matrix Three dimensional problems often require the solution of a 3 x 3 eigenproblem to generate principal values and directions FEAP includes a special routine to calculate the values and vectors for symmetric arrays The routine is used by a call to the subprogram as CALL EIG3 V D ROT CHAPTER 6 UTILITY ROUTINES 71 On call to the routine V 3 3 is a REAL 8 array containing the symmetric array to be diagonalized On return the eigenvalues are contained in D 3 and the vectors for each value in the columns of the V array A Jacobi method is used with ROT an integer parameter returning the number of rotations to diagonalize The routine is quite efficient compared to any attempt to compute vectors after closed form solution of the cubic for roots In addition to the general eigensolution above FEAP includes options to compute principal values of a symmetric second order tensor for two and three dimensional problems In two dimensional use the call to CALL PSTR2D SIG PV is used where SIG 4 stores stresses in the order 011 022 033 012 and returns principal values and directions in PV 3 in the order or 02 and 0 where the angle is in degrees between x and the 1 axis This routine does not use SIG 3 In three dimensions the principal v
63. oints and SW 5 is an array which returns three area coordinates and the quadrature weight SW 1 1 returns the volume coordinate L as defined in 2 3 SW 2 1 returns the volume coordinate La 1 SW 3 1 returns the volume coordinate L31 SW 4 1 returns the volume coordinate SW 5 1 returns the weight Wy Table 6 2 describes the admissible types number and location of quadrature points 6 2 Shape function subprograms Finite element approximations commonly use shape function subprograms to perform computations of the functions and their derivatives at preselected points often the quadrature points FEAP includes options to obtain the shape functions for some low order elements linear and quadratic order in one and two dimensions and linear shape functions for three dimensions In addition a cubic Hermitian interpolation routine is available The calling arguments for routines is summarized below CHAPTER 6 UTILITY ROUTINES 67 6 2 1 Shape functions in one dimension Lagrangian interpolation in one dimensional isoparametric forms may be obtained us ing the call CALL SHP1D S XL SHP NDM NEL XJAC where Parameter Description S Natural coordinate XL NDM Nodal coordinates for element NDM Spatial dimension of mesh NEL Number element nodes 2 or 3 SHP 2 NEL Shape function and derivative XJAC Jacobian transformation The shape functions are evaluated as SHP 1 i shape function deriv
64. ome large thus reducing the size of problem that FEAP can solve 5 6 2 Accessing history data for each element As noted above the data for each element is contained in arrays whose first word is located at hr nh1 hr nh2 where nh1 and nh2 are pointers for tn t 41 respectively and at hr nh3 for that not associated with time note that there are values for each CHAPTER 5 ADDING ELEMENTS 50 only if non zero values are assigned to nh1 and or nh3 during the isw 1 task Any other allocated data follows immediately after each first word It is a users responsibility to manage what is retained in each variable type however the order of placing the t and t 41 data into the NH1 and NH2 arrays should be identical There are no provisions to store integer history variables separately from double precision quantities It is necessary to cast the integer data as double precision and move to the history location For example using the statement hr nh3 5 dble ivarbl saves the value for the integer variable ivarbl in the sixth word of the NH3 element history array At a subsequent iteration for this element the value of the integer would be recovered as ivarbl int hr nh3 5 While this wastes storage for integer variables experience indicates there is little need to save many integer quantities and thus it was not deemed necessary to provide for integer history variables separately Although users may define new values for a
65. on has the command line GETData VALUes 35 is developed in the user function UMESH1 the first argument GETData is assigned to the variable uct as the name of the function see above however by default the other two words are ignored To recover their values the statement include chdata h CHAPTER 4 USER FUNCTIONS 17 character lct 15 real 8 ctl integer itl is added to the user function and the items recovered in the else option of the function using the statements lct yyy 16 30 call setval yyy 31 45 15 ctl If FEAP parameters or expressions are used to define a numerical entry the subprogram setval will make the appropriate evaluation If necessary the real value for ctl can be cast into an integer using itl nint ctl Note that the total number of added words must be 15 items or less this is imposed by the total of 16 items on any FEAP input record 4 2 User Material Models Users may add material models to elements by appending subprograms UMATIn and UMATLn where n have values from 0 to 9 to the FEAP system The subprogram UMATIn defines parameters used by the model and the subprogram UMATLn is called by the element for each computation point i e the quadrature point receives the value of a deformation measure as input and must return the value of stress and tangent moduli as output To activate a user material model the input data for the mesh MATErial command must include a statement with UCON
66. parameters It is highly recommended that users use include files rather than giving equivalent common statements directly If later releases of the FEAP program revise contents in a common block it will only be necessary to recompile the user routine rather than change all the common statement definitions Also by defining the correct path in the makefile in or your compiler it is not necessary to modify any routine when switch from 32 bit machines to 64 bit machines Chapter 2 DATA INPUT AND OUTPUT FEAP includes utilities to perform input and to output small arrays of data Users are strongly encouraged to use the input utilities but often may wish to use their own utilities to output data 2 1 Parameters and Expressions The subroutines PINPUT and TINPUT are input subprograms used by FEAP to input each data record They permit the data to be in a free form format with up to 255 characters on each record as well as to employ expressions parameters and numerical representations for each data item These routines also should be used to input data in any new program module developed The PINPUT routine returns data to the calling subprogram in a double precision array The following statements may be included as part of the routine performing the input subroutine xxx logical errck pinput integer ior iow ilg common iofile ior iow ilg real 8 td 5 1 if ior lt 0 write 3000 errck pinput td 5 if errck go to 1
67. pppcol color switch where color is an integer defining the color number and switch should be zero The color values are given in Table 4 1 Chapter 5 ADDING ELEMENTS FEAP permits users to add their own element modules to the program by writing a single subprogram called subroutine elmtnn d ul xl ix tl s r ndf ndm nst isw where nn may have values between 01 and 50 Each element subprogram must be added before loading the FEAP library since dummy subprograms are included in the library to avoid unsatisfied externals The basic structure for an element routine is shown in Figures 5 1 and 5 1 Information is provided to the element subprogram through data passed as arguments and data passed in common blocks The data passed as arguments consists of eleven subroutine elmtnn d ul xl ix tl s r ndf ndm nst isw c Prototype FEAP Element Routine nn 01 to 50 implicit none c Common blocks See Figure 5 2 integer ndf ndm nst isw integer ix nen1 1 real 8 d ul ndf xl ndm tl s nst nst r nst if isw eq 0 and ior lt 0 then c Output element description write Elmt 1 Element description Figure 5 1 FEAP Element Subprogram Part 1 27 CHAPTER 5 ADDING ELEMENTS 28 elseif isw eq 1 then Input output of property data after command mate d stores information for each material set Return nhi number of nh1 nh2 words element Return nh3 number of nh3 words element 0000 elseif isw eq 2 then C
68. r node ndm Space dimension of mesh nst Size of element arrays S and R N B Normally nst ndf nen isw Task parameter to control computation See prototype element in Figure 5 1 Table 5 1 Arguments of FEAP Element Subprogram CHAPTER 5 ADDING ELEMENTS 30 It is not necessary to implement all other tasks in an element however for those tasks that are not implemented it is important that the element routine not perform any calculations Thus if the form of the branch is programmed as an IF THEN ELSE construct as shown in Fig 5 1 then the ELSE should not carry out any operations unless all options for ISW are programmed Similarly if the element is programmed using a SELECT CASE form as isw Task Type Description Access Command 0 O Output label SHOW ELEM 1 R Input d parameters Mesh MATE n 2 O Check elements Soln CHECk 3 R Compute tangent residual Soln TANG Store in S r UTAN 4 O Output element variables Soln STRE 5 E Compute cons lump mass Soln MASS Store in S r MASS LUMP 6 R Compute residual Soln FORM REAC Plot REAC 7 Surface load tangents Mesh SLOAd 8 O Nodal projections Soln STRE NODE Plot STRE PSTR 9 O Damping Soln DAMP 10 O Augmented Lagrangian update Soln AUGM 11 O Error estimator Soln ERRO 12 R History update Soln TIME For special history treatments else return 13 O Energy momentum Soln TPLO ENER 14 R Initialize history BATCh INTEr 15 O Body force Mesh BODY
69. r users to prepare plot outputs unrelated to elements The plot command PLOT UPLOt vi v2 v3 initiates a call to the subroutine UPLOT which has the basic structure SUBROUTINE UPLOT CT IMPLICIT NONE REAL 8 CT 3 END The argument CT contains the values for the three parameters v1 v2 v3 The default color is white Direct plots in screen coordinates lower left at 0 0 upper right at 1 1 may be given using the statement CALL DPLOT XS YS IP where XS YS are between zero 0 and one 1 and IP is interpreted according to Table 6 4 Panels are closed using CALL CLPAN and colors treated according to values specified in calls to PPPCOL Bibliography 1 G G Weber A M Lush A Zavaliangos and L Anand An objective time integration procedure for isotropic rate independent and rate dependent elastic plastic constitutive equations International Journal of Plasticity 6 701 744 1990 21 O C Zienkiewicz and R L Taylor The Finite Element Method volume 1 McGraw Hill London 4th edition 1989 3 O C Zienkiewicz and R L Taylor The Finite Element Method The Basis vol ume 1 Butterworth Heinemann Oxford 5th edition 2000 4 M Abramowitz and I A Stegun editors Handbook of Mathematical Functions Dover Publications New York 1965 5 T J R Hughes The Finite Element Method Linear Static and Dynamic Analysis Prentice Hall Englewood Cliffs N J 1987 6 M Gellert and R Harbord Moderate degree c
70. s pata ES pode e dn ae a 36 Material Parameters cyst des e Bek a at 37 Tangent Parameters ara ba abs es ee s es s 41 Element Plot Definition Subprograms 47 Momenta and Energy Assignments 54 Quadrature for triangles code fare Sie gs hb Bod ANN BA e ge abe 65 Quadrature for tetrahedra 2 la is rs D S R 66 Color pallet fer FEAP plots war pi maarssen ween ea eh 74 Values for control of plots ar 20 3 0 A ELA Ake Moa nag ak 75 Chapter 1 INTRODUCTION In this part of the FEAP manual some of the options to extend the capabilities of the program are described We begin by describing the utilities provided in FEAP for use in data input Options to add user commands for mesh and command language extensions is then described and finally the method to add an element to the program is described 1 1 Setting Program Options The size of problems which may be solved by FEAP depends on the amount of memory available in the computer as well as solution options used Memory for the main arrays used to solve problems is dynamically allocated during the solution Arrays are allocated and deallocated using a system subprogram PALLOC or for user developed modules using subprogram UALLOC Further information on use of these routines is given in Section 3 The IPR parameter in the FEAP82 module controls the specification of the size of REAL variables For typical UNIX and PC systems all real variables s
71. s which are ob tained from the negative of the terms multiplying the nodal displacements Accord ingly Ral _ 1 0 5 1 Az 1 4 Uj naar eead Al 6 e ds Hal is the residual for the i coordinate direction For constant properties and loading over an element length note that for this case the stress will also be constant since strains are constant on the element the above may be integrated to yield CHAPTER 5 ADDING ELEMENTS 45 Ar Ri 1 1 Oss pAL 2 1 For the present we assume the material model is a linear elastic in which the stress is related to strain through Oss E Ess where E is the Young s modulus Based on a linear elastic material the term in the residual involving o may be written Ge SRE BA aaa us L 12 g lsin St Ba j 1 32 as For the linear elastic material a stiffness matrix may be expressed as sa EA Ax H kij Shy where HA kij Ts Az 2 The residual may now be written using a stiffness and mass matrix as d m Ra 1 1 kij kij Mii Mia s Bd 277 A 3 F kij ma Maa ibi ee with pAL pAL Mii Maa gt 5 Th Mo lt lt 3 6 For non linear material behavior the residual must be computed using Equation 5 1 with the stress replaced by the value computed from the constitutive equation The integration method for nodal quantities is taken as Newmark s method described in Section 5 2 The resid
72. sed in size only when the new size is less than the old size by the assigned value The last parameter which may be set in the FEAP82 module is the level for displaying available commands when the HELP command is used while in mesh solution or plot CHAPTER 1 INTRODUCTION 3 mode FEAP contains a large number of commands which are not commonly used by many users To control the default number of commands displayed to users the com mands have been separated into four levels 0 Basic 1 Intermediate 2 Advanced and 3 Expert The level to be displayed when using the HELP command is given may be set in the integer variable HLPLEV That is setting hlplev 1 Intermediate results in commands up to the intermediate level being displayed It is possible to raise or lower the level during execution using the command MANUal level where level is the numerical value desired When developing program modules it is often desirable to have output of specific quantities available e g tracking the change in some parameters during successive iterations FEAP provides for a switch to make the outputs active or inactive during an execution The switch is named debug and placed in integer ndebug logical debug common debugs ndebug debug The value of the debug is set true by the solution command DEBUg and false by the command DEBUg OFF Thus placing code fragments into modules as if debug then write iow LABEL list wr
73. set to NH1 2 history variables default 0 Offset to NH3 history variables default 0 Number history variables element NH3 Finite rotation update number for UROTxx Get tangent from element if 0 if gt 0 numerically differ entiate residual to obtain tangent Equation number for element Lagrange multiplier Table 3 5 Element control array IE use for material number ma 12 CHAPTER 3 ALLOCATING ARRAYS 13 The subprograms PALLOC and UALLOC may also be used to destroy a previously defined array This is achieved when the length of the array is specified as zero 0 For example to destroy the arrays defined as TEMP1 and TEMP2 the statements setvar palloc 111 TEMP1 0 2 setvar 1 palloc 112 TEMP2 0 are given Use of these statements results in the pointers np 111 and np 112 being set to zero and the space used by the arrays being released for use by other allocations at a later point in the program A call to PALLOC or UALLOC for any previously defined array but with a different non zero length causes the size of the array to be either increased or decreased For user defined arrays specified in UALLOC care should be exercised in selecting the alphanumeric NAME parameter which is limited to 5 characters so that conflicts are not created with existing names use of the SHOW DICT command is one way to investigate names of arrays used in an analysis or check the names already containe
74. turned if 2 quadratic interpolation if 3 the interpolation is generated plus a cubic bubble in the seventh function Giving the IORD parameter as a negative returns hierarchical form for mid side nodes 6 2 3 Shape functions in three dimensions Three dimensional Co isoparametric interpolation on bricks of linear order i e 8 node elements may be obtained using the subprogram call CALL SHP3D SS XJAC SHP XL NDM NEL where CHAPTER 6 UTILITY ROUTINES 70 Parameter Description SS 3 Natural coordinates n XL NDM 8 Element coordinates in local order NDM Spatial dimension mesh 2 or 3 NEL Number nodes on element 8 linear brick 20 serendipity quadratic 27 Lagrangian quadratic 64 Lagrangian cubic SHP 4 8 Shape functions and derivatives XJAC Jacobian transformation from xyz to ENC The array SHP stores the values in the order SHP 1 i derivative with respect to x SHP 2 i derivative with respect to y SHP 3 i derivative with respect to z SHP 4 i shape function Three dimensional Co isoparametric interpolation on tetrahedra of linear order i e 4 node elements may be obtained using the subprogram call CALL TETSHP SS XL NDM NEL XJAC SHP where Parameter Description SS 4 Volume coordinates Li La Ls La XL NDM 4 Element coordinates in local order NDM Spatial dimension mesh 3 NEL Number of nodes on element 4 10 11 14 15 XJAC Jac
75. ual and tangent matrix for a Newton type method are now available and may be inserted into R and S after noting that for the truss that the damping matrix C is zero The residual may be programmed directly from Equation 5 1 and an implementation using the two dimensional form r ndf nen is shown in Figure 5 5 Similarly using the results from Section 5 2 the tangent matrix for the truss may be programmed as indicated in Figures 5 6 and 5 7 CHAPTER 5 ADDING ELEMENTS 46 5 4 Additional Options in Elements FEAP permits some additional options to be included within element tasks 5 4 1 Task 1 Options Often it is necessary to use several element types to perform an analysis For example it may be necessary to use both truss and frame bending resistant elements to perform an analysis As developed in Section 5 3 the truss element has one degree of freedom for each spatial dimension whereas the frame element must have additional unknowns to represent the bending behavior For nodes connected only to truss elements it is not necessary to have the additional degrees of freedom active and a user would be required to specify restraint conditions for these nodes and degrees of freedom By inserting the following lines of code into the truss element routine for the isw 1 task FEAP will automatically eliminate any unneeded degrees of freedom do i ndmti ndf ix i end do i Note that for isw 1 the ix parameter is not used to pass th
76. ubature formulas for 3 D tetrahedral finite element approximations Communications in Applied Numerical Methods 7 487 495 1991 77 Index Adding elements 27 ELMTnn subprogram 27 Arguments on subprogram ELMTnn 29 Array allocation 8 PALLOC subprogram 9 Removing arrays 10 UALLOC subprogram 9 Color in plots 26 Common block definitions 33 Common statements 3 Data input Converting text to real 7 PINPUT 5 TINPUT 6 Data output Integer arrays 7 Real arrays 7 Debug output 3 Eigenvalues for 3 x 3 matrix 70 Element array names 11 Element common block variables 57 Element connection array 12 Element control array 12 Element energy computation 53 Element history variables 48 Accessing for each element 49 Element options Projection to nodes 47 Task 1 ISW 1 46 Task 6 ISW 6 46 78 TPLOt definitions 46 Element plots User 73 Element task definitions 30 Finite rotation 50 Help files 1 Include statements 3 Locating arrays PGETD 13 Material models 17 Auto time step control 21 History variables 20 Internal variables 20 UCON use 17 UMATIn subprogram 18 UMATLn subprogram 20 Material property variables 31 Mesh array names 11 Mesh input UMESHn subprogram 15 Mesh manipulation UMANI subprogram 23 Mesh plots 71 Brick element 73 Line element 72 Point element 72 Quadrilateral element 73 Tetrahedral element 73 Triangular element 72 INDEX Newto
77. ure of the algorithms employed by FEAP to solve problems Each problem to be solved using an ELMTnn routine is established in a standard finite element form as described in standard references e g The Finite Element Method 4th ed by O C Zienkiewicz and R L Taylor McGraw Hill London 1989 vol 1 1991 vol 2 Here it is assumed this step leads to a set of non linear ordinary differential equations expressed in terms of nodal displacements velocities and accelerations given by u t t and u t respectively We denote the differential equation for node i as the CHAPTER 5 ADDING ELEMENTS 33 Variable Definition nh1 Pointer to t history data nh2 Pointer to t 4 history data nh3 Pointer to element history ior Current input logical unit iow Current output logical unit nph Pointer to global projection arrays ner Pointer to global error indicator erav Element error value j int J integral values ndf Maximum dof node ndm Mesh space dimension nent Dimension 1 on IX array nst Size of element matrix nneq Total dof in problem ttim Current time dt Current time increment ci Integration parameters hr Real array data mr Integer array data Table 5 4 FEAP common block definitions residual equation To solve for the nodal displacements velocities and accelerations it is necessary to introduce an algorithm to integrate the nodal quantities in time specify a constitutive relation
78. used to perform trans formations or manipulations on previously prescribed data These commands appear between the mesh input END command and the first INTEractive or BATCh solution com mand To add a mesh manipulation command a subprogram with the name UMANIn where n has a value between 0 and 9 must be written compiled and linked with the program The basic structure of the routine UMANT1 is subroutine umanil prtu c User defined routine to manipulate mesh data for FEAP implicit none include umaci h 1 Contains UCT variable logical prtu Set name 1 to user defined if pcomp uct mani 4 then uct Set user defined command name else User execution function statements follow end if end The parameter PRTU is a logical parameter which is set to false when the NOPRint mesh command is given and to true when the PRINt command is used default is true The common block UMAC1 transfers the character variable UCT for the name of CHAPTER 4 USER FUNCTIONS 24 the command The default names are MANn where n is the same as the routine name number Assignment of a unique character name which must not conflict with names already assigned for mesh input commands should be used to replace the xxxx shown After FEAP completes the input of mesh data it scans all of the UMANIn routines and replaces the command names 1 etc by the user furnished names The ability to get array names as shown in C
79. uted as p f ova Q n l x x pd Q the angular momenta as the kinetic energy K p v v dQ Q J W C 40 ye E I Fa dr the stored energy as U and the work by external loads as The value of the displacement and velocity at the current time t are passed in ul i j 1 and ul i j 4 respectively Note that this is true no matter which time integration algorithm is specified 5 9 A Non linear Theory for a Truss A simple non linear theory for a two or three dimensional truss which may undergo large displacements for which the strains remain small may be developed by defining the axial strain approximation in each member as CHAPTER 5 ADDING ELEMENTS 59 where un is a displacement component normal to the axis of the member The virtual strain from a linearization of the strain is given as ats 00U OUnj eso F j 1 An algorithm to define the two orthogonal unit vectors which are normal to the member may be constructed by taking v EL where k is a direction for which a minimum value of the direction cosine l exists for a 2 dimensional problem defined in the z z plane v may be taken as ez Now v xl n Tv xd and n lxn Using these vectors the two normal components of the displacement are given by d uu O nj ufs t Y nauds t i l and the derivative by OUn 2 Ou Collecting terms and combining with previously defined qua

FEAP - - A Finite Element Analysis Program

Contents

Download Pdf Manuals

Related Search

Related Contents