Home

FEAP - - A Finite Element Analysis Program

1. 32 5 4 Node Truss Element 2 4 094k oh KABA aS SO 43 5 5 Element residual for two node truss 46 5 6 Truss Tangent Matrix Parts 3 oe a ES ace SS ARES 47 5 7 Truss Tangent Matrix Part aa a daa Tone Ang 48 5 8 Element variable projection routine 52 111 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 00 11 Solution Array Names Numbers and Sized 0 11 Element Array Names Numbers and Sizes 11 Element connection array IX use forelemente 12 Element control array IE use for material number ma 12 Color Table for Plots os dos an GA ae s k sas usa Ma es ed oe old es 25 Arguments of FEAP Element Subprogram 28 Task Options for FEAP Element Subprogram 29 FEAP common block definitions 2 0 000 33 FEAP common block definitions 2 5 2 4 wf paka ha ka KG WG 34 Material Parameters 1 2 ite AN Ale At a 35 Material Parameters sana yam lana lee Se Are SoS 36 Material Parameters 202 0 es 2 ble ee BO GD LA na gamo IN LEE s 37 Material Parameters Ana plas pata Be he es es Nga LG ei 38 Tangent Parameter e s aus Yow ee kna BE NAAN DAWA 41 Element Plot Definition Subprograms 50 Momenta and Energy Assignments a 55 Quadrature for triangles LIO AAA AA
2. 7 XL NDM NEL Element coordinates in local order NDM Spatial dimension mesh 2 or 3 NEL Largest local node number on element 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 z y to n The array SHP stores the values in the order SHP 1 i derivative with respect to or x SHP 2 1 derivative with respect to y or y SHP 3 i shape function Two dimensional Co isoparametric interpolation on triangles of linear or quadratic order may be obtained using the subprogram call CALL TRISHP SS XL NDM IORD XJAC SHP where Parameter Description 553 Area coordinates L1 La Ls XL NDM x Element coordinates in local order NDM Spatial dimension mesh 2 or 3 IORD Order of interpolation 1 2 or 3 XJAC Jacobian transformation from z 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 y 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 returned 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 no
3. which is placed in the file comblk 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 a 4 GByte range The pointers in the integer8 subdirectory have 64 bit lengths and should be used for larger 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 When using 64 bit pointers users must be careful to always define the address of an array in a calling statement to also be 64 bits in length For example use of CHAPTER 3 ALLOCATING ARRAYS 9 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 integers Using this scheme permits direct reference to either real 8 or integer a
4. 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 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 4 USER FUNCTIONS Number Color CONDO gt ON HO O 10 11 12 13 14 15 16 17 18 Black White Red Green Blue Yellow Cyan Magenta Orange Coral Green Yellow Wheat Royal Blue Purple Aquamarine Violet Red Dark Slate Blue Gray Light Gray Table 4 1 Color Table for Plots 25 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 elm
5. as Fung model energy parameter 175 as Fung model energy parameter 176 ag Fung model energy parameter 177 a7 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 Mg mass scaling factor 184 c Estimate on maximum wave speed 185 Augmentation switch 186 200 Unused Table 5 8 Material Parameters CHAPTER 5 ADDING ELEMENTS 39 example consider the case where the solution is parameterized in terms of increments of the displacements i e a is the displacement vector u For this case the tangent matrix is we do not show dependence on the iteration k for simplicity of notation S du _ ORG iu _ OR OV _ OR Jax Uu Bea hg Ad Ov Ou da Ou 7 Note that from the Newmark formulas A s OVEN NOUR OBE a A Ou BAR du da du BAt in which 94 is the Kronnecker delta identity matrix for the k j nodal pair From the residual we observe that K ae ga OR M OR da du Cd Ov define the tangent stiffness damping and mass respectively Thus for the Newmark algorithm the total tangent matrix in terms of the incremental displacements is 1 M BAP 4 po T Sy Kg BALCH For other choices of increments the tangent may be written in the general form S GK amp y c3Miy where the c are scalar quantities involvin
6. ci Integration parameters hr Real array data mr Integer array data Table 5 4 FEAP common block definitions 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 define the relationships between the increments Note that only scalar multipliers involving 0 y and At are involved between the different measures The tangent matrix for the transient problem using Newmark s method may be ex pressed in terms of the incremental displacement velocity or acceleration As an CHAPTER 5 ADDING ELEMENTS Parameter Name Description 1 E Young s modulus DI u Poisson ratio 3la Thermal expansion coefficient 4 p Mass density 5 Quadrature order for arrays 6 Quadra
7. 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 recoverd 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 function 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 character lct 1i5 realx8 ctl integer itl is added to the user function and the items recovered in the else option of the function using the statments CHAPTER 4 USER FUNCTIONS 16 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
8. returns the weight W Table 6 2 describes the admissible types number and location of quadrature points CHAPTER 6 UTILITY ROUTINES 61 Type Number Location Points 1 Centroid O h 4 Interior O h 11 Interior O h4 16 Interior O h He WN e Table 6 2 Quadrature for tetrahedra 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 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 derivative along the axis of the element and SHP 2 i the shape function N In calculations integrals are r
9. where the entries in each submatrix are given as Ry Ea Rq The two dimensional form places the entries R as columns Accordingly R Ri R R3 el The two forms for defining the residual r are eguivalent based on the Fortran ordering of information into double subscript arrays If ndf is larger than needed for the element 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 xl i j 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 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 CHAPTER 5 ADDING ELEMENTS Al Parameter Description ctan 1 C1 Multiplier of s matrix for ul i j 1 terms e g stiffness matrix multiplier ctan 2 c2 Multiplier of s matrix for ul i j 4 terms e g damping matrix multiplier ctan 3 c3 Multiplier of s matrix for ul i j 5 terms e g mass matrix multiplier Table 5 9 Tangent Parameters to add these
10. 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 EA eps Element axial force tt 2 e ps Element axial strain Y ba 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 1 saves force for element nn stress nn 2 saves strain for element mn 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 values Accordingly a Ng Ga CHAPTER 5 ADDING ELEMENTS 51 where g is any value which is to be projected to nodes e g a stress or strain N are shape functions for the element type considered and o 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 m
11. 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 subprogram 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 5 2 Non linear Transient Solution Forms Before describing the steps in developing an element we summarize first the basic structure 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 equat
12. Auto time step control 19 History variables 19 Internal variables 19 UCON use 16 UMATIn subprogram 17 UMATLn subprogram 17 Material property variables 30 Mesh array names 11 Mesh input UMESHn subprogram 14 Mesh manipulation UMANIn subprogram 21 Mesh plots 65 Brick element 66 Line element 66 Point element 66 Quadrilateral element 66 Tetrahedral element 66 Triangular element 66 Newton solution Element residual array storage 40 Element tangent array storage 39 Non linear problems 32 Tangent array definition 34 Numerical integration 58 One dimension 58 Three dimension 60 Two dimension 59 INDEX Plot command UPLOT subprogram 69 UPLOTn subprogram 23 Plot utility function Fills 67 Lines 67 Plot utility functions Color 24 Fills 24 Lines 24 Precision Integer 1 Real 1 Principal stress Three dimensional 65 Two dimensional 65 Problem size 1 Quadrature 58 One dimension 58 Three dimension 60 Two dimension 59 Setting colors 24 Setting options Help level 2 Parsing statements 2 Plot prompts 2 PostScript 2 Shape functions 61 Brick isoparametric 64 One dimensional cubic Hermitian 62 One dimensional isoparametric 61 One dimensional natural derivative 62 Quadrilateral isoparametric 63 Tetrahedral isoparametric 64 Three dimensional isoparametric 64 Trianglular isoparametric 63 Two dimensional isoparametric 63 Solution array names 11 Solution command UMACRn subp
13. Az Az be i 1 Using the above displacement interpolations the variational equation for the truss may be expressed in matrix form as 011 un u z 4 5 Fe os Ads f i Pla pa 1 2G ds EN J ea bs bth le FEAP constructs the finite element arrays from the element residuals which are ob tained from the negative of the terms multiplying the nodal displacements Accord ingly po Bal 1 l A E bas 1 Na 1 4 Uj 5 12 ouad Jean g as i CHAPTER 5 ADDING ELEMENTS 45 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 _ Ra lt A 1 Oss Ax pAL 2 1 Ua R al 582 o OL E 6 1 21 a c 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 ds may be written as L Ati J Ati D ai ai P j 1 1 For the linear elastic material a stiffness matrix may be expressed as _ FA Az _ ky Shy oo Ee aa kaa pa nad where EA kij ES Az Az The residual may now be written using a stiffness and mass matrix as d ae Ra x 1 1 kij kij Uji Mai Mia s Ha Bd x gt B gt k kij Uj2 ma M 2 ibi
14. 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 ds Plane stress recovery 91 d3 Plane stress recovery 92 a3 Plane stress recovery Table 5 7 Material Parameters 37 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 ng Reference vector parameter 100 Cross section shape 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 Too 129 Reference absolute temperature 130 nseg Number of hardening segments 131 148 Segment data sets pYisoH kin 149 Total variables on frame section 150 Piezoelectric flag 151 159 Piezoelectric data 160 Initial stress flag 161 166 o Initial stresses constant 167 Tension compression only indicator 170 C Fung pseudo elastic model modulus 171 aj Fung model energy parameter 172 a2 Fung model energy parameter 173 a3 Fung model energy parameter 174
15. default ma Offset 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 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 contained 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
16. 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 i e values are reported for time tn 1 unless generalized midpoint integration methods are used For the present we will assume a is 1 CHAPTER 5 ADDING ELEMENTS 28 FEAP carries out tasks according to the value of the parameter ISW passed as the 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 0 to 6 8 and 10 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 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 It is not necessar
17. ee ee a 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 kr PAKA GA a a a E implicit none include iofile h logical pcomp character type 15 integer ni n3 real 8 vv 5 ud c Specify type of user model if pcomp type mati 4 then type E 1d Specify new name for model c Input output user data and save in ud array c Set values of n1 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 19 the trace of the strain in the parameter theta Thus 0 Eii 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 quad
18. if pcomp uct mac0 4 then uct xxxx else c User command statements are placed here endif end 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 uplot0 ct1 uprt c User plot command function implicit none include umaci h Contains the variable UCT logical uprt real 8 ct1 3 Set command word if pcomp uct p1t0 4 then uct xxxx 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 CHAPTER 4 USER FUNCTIONS 24 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 call plotl x1 x2 x3 ipen where x1 x2 x3 are coordinates values of the mesh The value of ipen ranges from
19. n 8 compro 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 Lisa to IX nen e Global node nen IX nent1 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 Nay 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 TE nie 8 ma Number history variables element NH1 and NH2 Element material type number ELMTO1 1 etc Element material type identifier
20. 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 Integer assignment float td 2 Realx4 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 to 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 is equivalent to errck pinput td nn CHAPTER 2 DATA INPUT AND OUTPUT 7 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 out
21. 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 Atnew 1 25 At 0 5 lt rmeas lt 0 8 Atnew At rmeas 0 8 lt rmeas 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 At 3 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 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 umani prtu c User defined routine to manipulate mesh data for FEAP implicit none include umaci h Contains UCT variable logical prtu CHAPTER 4 USER FUNCTIONS 22 Set name mani to user defined if pcomp uct mani 4 then uct xxxx Set user defined command name else c User execution function statements follow end if end The paramet
22. the current material number Care should be taken to minimize the number of history variables since for very large problems the memory requirements can become 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 tn 1 respectively and at hr nh3 for that not associated with time note that there are values for each 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 1 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 CHAPTER 5 ADDING ELEMENTS 54 While this wastes storage for integer variables experience indicates there is little need to save man
23. 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 Cior 1t 0 write 3000 errck pinput td 5 if errck go to 1 The parameters defined in the common block are ior input file unit number if negative input CHAPTER 2 DATA INPUT AND OUTPUT 6 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
24. 2 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 By default it is assumed that all plots are for two dimensional elements of quadrilateral shape and with between 4 and 9 nodes numbered as shown in the manuals vertex nodes first mid side nodes center node If other shapes are to be represented it is necessary to CHAPTER 6 UTILITY ROUTINES 66 include drawing instructions in the element routine segment where material properties are input i e ISW 1 A parameter describing the type of element and named TEL is contained in the labeled common ELDATA and thus this common must be included to activate the plot option call 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 3 Triangular element with 3 nodes obtained by call CALL PLTRI3 IEL 4 Triangular element with 6 nodes obtained by call CALL PLTRI6 IEL 5 Quadrilateral element with 4 nodes obtained by call CALL PLQUD4 IEL 6 Tetrahedral element with 4 nodes obtained by call CALL PLTET4 IEL 7 Brick element with 8 nodes obtained by call CALL PLBRK8 IEL Using these and internal extraction of element surfaces t
25. A AS e A e 60 Quadrature for tetrahedra 2 2 ab kaon oy Ken a e e r 61 Color pallet for FEAP plots oia o e aa de 68 Values for control of plots al Sal aa dei ol bene o 68 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 FEAP81 module controls the specification of the size of REAL variables For typical UNIX and PC systems all real variables should 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 paramete
26. APTER 5 ADDING ELEMENTS if isw eq 3 then c Compute element length L2 0 0d0 do i 1 ndm L2 L2 xl i 2 xl i 1 xx2 end do L sgrt L2 c Form stiffness multiplier dd ctan 1 xEAxL c Compute strain displacement matrix Lr 1 d0 L2 do i 1 ndm bb i 1 x1 1 2 x1 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 47 CHAPTER 5 ADDING ELEMENTS 48 c Compute stiffness terms N B ndm lt or ndf il O do ii 1 2 jl 0 do jj 1 2 do i 1 ndm do j 1 ndm s itit j j1 db i ii bb j jj end do end do ji jl ndf end do il 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 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 49 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 re
27. CON record and these CHAPTER 4 USER FUNCTIONS 17 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 values as vi 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 ni Currently t
28. FEAP A Finite Element Analysis Program Version 8 1 Programmer Manual Robert L Taylor Department of Civil and Environmental Engineering University of California at Berkeley Berkeley California 94720 1710 E Mail rltGce berkeley edu January 2007 Contents Introduction 1 1 Setting Program Options __ _ _ SA eet SoA td eC 1 2 Uses of Common and Include Statements Data Input and Output 2 1 Parameters and Expressions 2 ne esd rd Scere OS 22 Array Outp ts O A NA ai suqus ck Allocating Arrays User Functions 41 Mesh Input Functions UMESHn 00 0000 aa 4 1 1 Command line data e 42 User Material Models O 421 The UMATIn Module 4 2 2 The UMATLn Module 4 2 3 Auto time step COMO 2 AN an ea 4 3 Mesh Manipulation Functions UMANIn 4 4 Solution Command Functions UMACRn 0 4 5 Plot Command Functions UPLOTn 0 0 0 Adding Elements 5 1 Material property storage a td as ee Geet 5 2 Non linear Transient Solution Forms 04 5 3 Example 2 Node Truss Element 000 5 Hk Theory lora Truss s Lal A e hed Se ha PING 5 4 Additional Options in Elements 5 4 1 Task 1 Options e e a AAH S42 Task a s IRR AK ak ee ee A ed ee 5 5 Projection of element variables to nodes 0 5 6 Elements with Histo
29. M 1 Input d parameters Mesh MATE n 2 Check elements Soln CHECk 3 Compute tangent residual Soln TANG Store in S r UTAN 4 Output element variables Soln STRE 5 Compute cons lump mass Soln MASS Store in S r MASS LUMP 6 Compute residual Soln FORM REAC Plot REAC 7 Surface load tangents Mesh SLOAd 8 Nodal projections Soln STRE NODE Plot STRE PSTR 9 Damping Soln DAMP 10 Augmented Lagrangian update Soln AUGM 11 Error estimator Soln ERRO 12 History update Soln TIME 13 Energy momentum Soln TPLO ENER 14 Initialize history BATCh INTEr 15 Body force Mesh BODY 16 J integrals Soln JINT 17 Set after activation Soln ACTI 18 Set after deactivation Soln DEAC 19 NOT AVAILABLE used in modal base BASE Uses isw 5 in element 20 Element plotting Plot PELM 23 Compute element loads only ARCL 25 Zienkiewicz Zhu projection Soln ZZHU 26 Used to compute mesh boundary Called by default Table 5 2 Task Options for FEAP Element Subprogram CHAPTER 5 ADDING ELEMENTS 30 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 one Some of the options for additional data passed through common blocks is shown in Figure 5 2 with each variable defined in Table
30. R 6 UTILITY ROUTINES 60 Type Number Location Points 1 1 Centroid O h2 3 3 Mid sides O h3 3 3 Interior O h3 4 4 Interior O h 6 6 Nodal O h3 7 7 Interior O h9 Table 6 1 Quadrature for triangles 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 7 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 LIN 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 points 7 and SW 3 1 contains values of the points G and SW 4 1 contains values of the weights W Three dimensional quadrature on tetrahedra may be performed using the subprograms call CALL TINTSD C L LINT SW where L is a type indicator LINT returns the number of points and SW 5 is an array which returns three area coordinates and the quadrature weight SW 1 1 returns the volume coordinate L4 as defined in 2 3 SW 2 1 returns the volume coordinate La j SW 3 1 returns the volume coordinate L3 SW 4 1 returns the volume coordinate L41 SW 5 1
31. a panel is drawn with standard line drawing commands starting with CHAPTER 6 UTILITY ROUTINES 68 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 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 1 and continuing with a sequence 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 Ti 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 yi 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 paramete
32. 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 xy 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 xpointt1 The use of pgetd can lead 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 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 s
33. ae ee aa implicit none include cdata h Contains nen include strnum h 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 pa p i 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 blank common similarly data for t 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 nhi nh2 nh3 CHAPTER 5 ADDING ELEMENTS 53 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
34. 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 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 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 U
35. 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 FEAP81 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 decreased 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 FEAP81 module is the level for displaying available commands when the HELP command is used while in mesh solution or plot mode FEAP contains a large number of commands which are
36. ass lumping For example Ma f Na dQ Q defines a row sum form of projection 8 Using the above results in the set of equa tions and a least square fit with the finite element values amp gives the equation set Maa Ia aap Nod 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 quantities 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 t 1 and variables whi
37. ay be used in defining an acceptable value for rmeas The algorithm coded monitors the solution during a standard iteration process set by for example CHAPTER 4 USER FUNCTIONS 20 000000000000 aaa subroutine umatl1 eps theta td d ud hn hi nh istrt sig dd isw D a a 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 cr 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 paa A a Aer aa PE es ese yO cis Ie regis ene a Sie implicit none integer umat nh isw i realx8 td realt8 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 21 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
38. 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 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 comb1k h which is defined as realx8 hr integer mr common hr 1 mr 1 CHAPTER 1 INTRODUCTION 4 The dummy arrays hr 1 and mr 1 serve to pass all dynamically allocated arrays between subprograms using a pointer array contained in the common np and located in the include file pointer h See Section 3 for more details on use of pointers All include files are located in the directory include 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 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
39. ch 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 be attempted Before calling the element routine for each element FEAP transfers the required history variable 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 3An implementation of the Zienkiewicz Zhu projection method is implemented using ISW 24 CHAPTER 5 ADDING ELEMENTS 52 c Cc 00000000 c subroutine slcn2d sig shp xsj sg lint nel nes p s a a oe TO e aa asna Purpose Project element variables to nodes Inputs sig nes Stresses at quadrature points shp nel Shape functions at quadrature points xsj Volume element at quadrature points sg 3 Gauss points 1 2 and weights 3 lint Number of quadrature points nel Number nodes on element nes Dimension of stress array Outputs p nen Weights for lumped projection s nen Integral of variables Se oe Sos a
40. des CHAPTER 6 UTILITY ROUTINES 64 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 Parameter Description SS 3 Natural coordinates 7 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 1 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 ORDER XJAC SHP where Parameter Description SS 4 Volume coordinates Ly La Ls La XL NDM 4 Element coordinates in local order NDM Spatial dimension mesh 2 or 3 ORDER Interpolation order 1 linear 2 quadratic XJAC Jacobian 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
41. e arrays IX and IE respectively 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 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 32 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
42. ed as p f ova Q m Cw x x pan Q the angular momenta as the kinetic energy K J p v v dQ Q CHAPTER 5 ADDING ELEMENTS 59 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 stored energy as i W C dQ V fx Eww an r and the work by external loads as The value of the displacement and velocity at the current time t i 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 8 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 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 Adu Odunj Ounj ba Nag os T 2 s 55 An algorithm to define the two orthogonal unit vectors which are normal to the member may be constructed by taking v EL CHAPTER 5 ADDING ELEMENTS 56 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 12 plane v may be taken as ez N
43. epresented as HN Ng de TING Nil XJACE d 6 5 CHAPTER 6 UTILITY ROUTINES 62 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 S Natural coordinate SHP 2 NEL Shape function and derivative NEL Number element nodes 2 or 3 where SHP 1 i contains N and SHP 2 i the shape function N Cubic Hermitian interpolation e g for use in straight linear beam elements given by w NP q NY wo N 0 N 0 6 6 is obtained using the call 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 w 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 Nix 2 SHPW 2 i SHPT 2 i are second derivatives e g Ni xz 3 SHPW 3 i SHPT 3 i are third derivatives e g Ni rx and 4 SHPW 4 i SHPT 4 i are shape functions e g N CHAPTER 6 UTILITY ROUTINES 63 6 2 2 Shape functions in two dimensions Two dimensional Cp isoparametric interpolation on quadrilaterals of linear or quadratic order may be obtained using the subprogram call CALL SHP2D SS XL SHP XJAC NDM NEL IX FLG where Parameter Description SS 2 Natural coordinates
44. er 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 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 man1 etc by the user furnished names The ability to get array names as shown in Chapter 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 umacr0 1ct ctl uprt c User solution command function implicit none include umac1 h Contains the variable UCT logical uprt character lct 15 real 8 ct1 3 CHAPTER 4 USER FUNCTIONS 23 Set command word
45. ere in addition to previously defined quantities we define CHAPTER 5 ADDING ELEMENTS 43 Figure 5 4 2 Node Truss Element d is the spatial dimension of the truss 1 2 or 3 e x are the Cartesian coordinates in the d directions L is the length of the truss member du is a virtual displacement in direction x e is an acceleration in direction z v e b is a loading in direction x per unit length and e des 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 1 where t is time u is the displacement vector with components us l is a unit vector along the axis of the member with direction cosines l defined by Ox _ vig Ti CHAPTER 5 ADDING ELEMENTS 44 and i Tj are the coordinates of nodes 1 and 2 respectively The displacement components are interpolated on the 2 node truss member as ui s t 1 ualt ult in which us W3 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 d Ou Ou Gaw a dl Os Os 3 1 Using the interpolations for the displacement components yields d 1 Ess T2 gt Ar Au i 1 where Ar Lig Li L L and Au Ug u Thus in matrix form the strain is d a Fa gt
46. evelop 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 Vi tn and A typical example for an integration algorithm for these discrete quantities is New mark s method where 1 and Viltn 1 Viltn AIA y a t Yak tau 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 implement a Newton method it is necessary to linearize the residual equation For FEAP the Newton equation may be written as R RD R SE G ga 0 i Oa 7 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 CHAPTER 5 ADDING ELEMENTS 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
47. g 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 In 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 ndf the maximum 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 Sip Si3 S Sa So So3 S31 S32 S33 where S is the sub matrix for nodal pairs i 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 CHAPTER 5 ADDING ELEMENTS 40 is passed to the subprogram as ndf Thus Si Sp Si 1 1 1 S Sx Sap S33 39 Sx S32 65 ij Sndf mdf in which se is an array coefficient for nodal pair z 7 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 R R2 R Rs
48. he parameter n3 may be set but is not available to the user material model Thus all history variables must be retained in the ni list Use of history variables is described later as part of the UMODEL module 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 P e11 22 E33 Y12 723 Y31 where yi 26jj 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 Current values of the strains are as mentioned above passed in the array eps 6 and CHAPTER 4 USER FUNCTIONS 18 subroutine umatil type vv ud n1 n3 Ce Pa aa i
49. he program is able to make some hidden surface plots in three dimensions Failure to include a plot sequence will usually result in unpredictable plots of the mesh and contour values CHAPTER 6 UTILITY ROUTINES 67 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 vi v2 v3 In interactive mode in the plot environment it is only necessary to enter 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 calls 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 A straight line segment may be drawn to the screen in the current color between the coordinates 1 y1 21 and x5 Y2 22 using the commands CALL PLOTL X1 Y1 Z1 3 CALL PLOTL X1 Y1 Z1 2 Here the basic command is 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
50. ions CHAPTER 5 ADDING ELEMENTS 31 character 4 o head common bdata o head 20 integer numnp numel nummat nen neq ipr common cdata numnp numel nummat nen neg ipr integer nstep niter naugm titer taugm tform common counts nstep niter naugm titer taugm tform integer iaugm iform intvc iautl nstepa common counts iaugm iform intvc iautl nstepa realx8 dm integer n ma mct iel nel common eldata dm n ma mct iel nel realx8 tt common elplot tt 200 realx8 bpr ctan psil common eltran bpr 3 ctan 3 psil real 8 ut common eluser ut 200 integer nhi nh2 nh3 nlm common hdata nhi nh2 nh3 nlm integer ior iow ilg common iofile ior iow ilg integer nph real 8 erav j int jshft common prstrs nph ner erav j int 3 jshft integer ndf ndm neni nst nneq ndl nnim nadd common sdata ndf ndm nen1 nst nneq ndl nnlm nadd realx8 ttim dt c1 c2 c3 c4 c5 chi common tdata ttim dt c1 c2 c3 c4 c5 chi realx8 hr integer mr common hr 1 mr 1000 Figure 5 2 FEAP Element Common Blocks expressed in terms of nodal displacements velocities and accelerations given by u t u t and t respectively We denote the differential equation for node i as the CHAPTER 5 ADDING ELEMENTS 32 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 and d
51. 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 33 where a is one of the variables at time t 4 e g Uj tn41 We define sw Bi jw K Oa and solve sda RP The solution is updated using k 1 _ k k a Aa da In the above k is the iteration number for the Newton algorithm To start the solution for each step FEAP sets AO bner a tn CHAPTER 5 ADDING ELEMENTS 34 Variable Definition nh1 Pointer to t history data nh2 Pointer to t 1 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
52. 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 CHAPTER 6 UTILITY ROUTINES 65 the values and vectors for symmetric arrays The routine is used by a call to the subprogram as CALL EIG3 V D ROT 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 c1 02 and 0 where the angle is in degrees between z and the 1 axis This routine does not use SIG 3 In three dimensions the principal values are obtained using the call CALL PSTR3D SIG PV where SIG 6 stores stresses in the order 041 022 033 012 023 031 and returns principal values in PV 3 in the order oi 0
53. n 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 TEMP1 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 the arrays using their pointers e g in a subroutine call as include pointer h include comblk h sara laca 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 th
54. nitial yield stress Mises 42 Yeo Final yield stress Mises 43 0 Exponential hardening rate 44 Hise Isotropic hardening modulus linear 45 Hkin Kinematic hardening modulus linear 46 Yield flag AT Bi Orthotropic thermal stress 48 Bo Orthotropic thermal stress 49 83 Orthotropic thermal stress 50 Error estimator parameter 51 yy Viscoelastic shear parameter 52 T Viscoelastic relaxation time 53 va Viscoelastic shear parameter 54 To 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 36 CHAPTER 5 ADDING ELEMENTS Parameter 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 71119 Ground acceleration factor 72 9 Ground acceleration factor 73 93 Ground acceleration factor 74 pi Ground acceleration proportional load number 75 po Ground acceleration proportional load number 76 p3 Ground acceleration proportional load number 77 ao Rayleigh damping mass ratio 78 ay Rayleigh damping stiffness ratio 79 Plate
55. 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 CHAPTER 1 INTRODUCTION 3 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 writes 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 which are used to pass parameters and small array values
56. ow vxl n v x1 and n x ni Using these vectors the two normal components of the displacement are given by d Unj 5 t nj u s t Y nguls t i 1 and the derivative by Ounj Ou ng Po la Collecting terms and combining with previously defined quantities the virtual strain may be written as du a El d 1 Sips j 1 ss where After differentiation of the displacement field the discrete form of the virtual strain is given by 1 O s3 q ou du A Substituting the above virtual strain expression into the weak form gives the modified residual expression ayal 1 Si L 2 1 a motos Ol aa at bib 6 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 TA g I 08 and Oou Odu A065 n On n gt n gt EA CHAPTER 5 ADDING ELEMENTS 57 If the n are constructed as column vectors then the tensor product becomes a matrix defined as T T G n Qn n n nini nn With these definitions the tangent matrix for the non linear problem is given as Ky Tl bo ol Notice that for the linear problem Az L gi thus the only difference between the linear and non linear problem is the definition of Ess in terms of displacements the modification for geometric effects for the g and the second term on the tangent matrix which i
57. plastic constitutive equations International Journal of Plasticity 6 701 744 1990 2 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 LA 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 70 Index Adding elements 26 ELMTnn subprogram 26 Arguments on subprogram ELMTnn 28 Array allocation 8 PALLOC subprogram 9 Removing arrays 10 UALLOC subprogram 9 Color in plots 24 Common block definitions 34 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 64 Element array names 11 Element common block variables 31 Element connection array 12 Element control array 12 Element energy computation 54 Element history variables 51 Accessing for each element 53 Element options Projection to nodes 50 Task 1 ISW 1 49 Task 6 ISW 6 49 TPLOt definitions 49 Element plots User 67 Element task definitions 29 Help files 1 Include statements 3 Locating arrays PGETD 13 Material models 16
58. present 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 0 end do i Note that for isw 1 the ix parameter is not used to pass the 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 50 Routine Name Description PLTLN2
59. put 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 parameters 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 num 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 realx8 hr integer mr common hr 1 mr 1
60. r i e variable CS is available in labeled common PVIEW In this case one can CHAPTER 6 UTILITY ROUTINES 69 add the statements assuming here that the displacements correspond to the coordinate directions DO NE 1 NEL DO I 1 NDM XP I NE XL I NE CS UL I NE END DO 1 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 for 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
61. r 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 this option to work properly it is necessary to define the path name to each manual CHAPTER 1 INTRODUCTION 2 page in the FEAP81 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 normally 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 FEAP81 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
62. rature 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 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 x 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 realx8 rmeas rvalu logical aratfl common elauto rmeas rvalu 3 aratfl The three parameters m
63. rogram 22 Table of colors 25 12 Transient solution CTAN definition 41 Newton method 32 Non linear 30 Truss element Newmark method 45 Non linear theory 55 Tangent modulus 56 Theory 41 Variational eguation 44 Weak form 44
64. rrays 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 array 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 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 CHAPTER 3 ALLOCATING ARRAYS 10 As an example for the use of the above allocatio
65. ry Variables 24 see GG Kamang ea AY 5 6 1 Assigning amount of storage for each element 5 6 2 Accessing history data for each element WwW al 14 14 15 16 17 17 19 21 22 23 CONTENTS li 5 7 Energy Computation na lt b ga gage A ra od pak dt acd 54 5 8 A Non linear Theory fora Truss 200 4 55 6 Utility routines 58 6 1 Numerical quadrature routines 2 58 6 1 1 One dimensional quadrature 58 6 1 2 Two dimensional quadrature 59 6 1 3 Three dimensional quadrature 60 6 2 Shape function subprograms Hang ah BA RES BAN 61 6 2 1 Shape functions in one dimension 61 6 2 2 Shape functions in two dimensions 63 6 2 3 Shape functions in three dimensions 64 6 3 Eigenvalues for 3 x 3 matrix ek SA Oe ee 64 Od PIGETOULINED are Di See SS Bee A a kta 65 6AT Mesh plots AA 65 6 4 2 Element data plots pel Slack ed Pech de 67 Oa Other ser plots iy oa Ka Lan a ade 69 List of Figures 4 1 Sample UMATIT module a a ak a ok 18 4 2 Sample UMATLn module Sa SB wR Keg ae ats 20 5 1 FEAP Element Subprogram PST aa ee ata 26 5 1 FEAP Element Subprogram Part 22 oa aang Bok BA LAWA 27 5 2 FEAP Element Common Blocks a 4 4 bon oi duta ee ee ee Be 31 5 3 FEAP Element Common Blocks using Includes
66. s 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 dimensional quadrature Two dimensional quadrature on quadrilateral domains may be performed by repeated one dimensional integration The two dimensional integrations are approximated by II Kemasan Y Kem wi 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 LIN 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 CALL TINT2D L LIN 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 Ls as defined in 2 3 SW 2 1 returns the area coordinate Ly SW 3 1 returns the area coordinate La SW 4 1 returns the weight W Table 6 1 describes the admissible types number and location of quadrature points CHAPTE
67. s sometimes called the geometric stiffness part 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 integration 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 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 Max F W 6 1 1 l 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 The Gauss Legendre formula has points which are all less than unity The subpro gram call CHAPTER 6 UTILITY ROUTINES 59 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 i
68. sw 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 CHAPTER 5 ADDING ELEMENTS 42 5 3 1 Theory for a Truss The governing equations for a typical truss member element shown in Figure 5 4 are the balance of momentum equation OA Ao Os the strain displacement equation for small deformations Ab pA _ Qu LOG and a constitutive equation For example considering a linear elastic material the constitutive equation may be written as Oss Bes 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 Es is a strain along the truss member axis and e Gss is the stress on a truss cross section The equations may also be deduced from the variational equation d d Ol doo Ads S swpainas S swtnds Oller L i 1 YL i 1 YL where dll contains the boundary and loading terms not associated with an element Wh
69. 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 returned 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 defined by i
70. tnn 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 A ul ndf x1 ndm t1 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 26 CHAPTER 5 ADDING ELEMENTS para 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 eg 2 then Check element for errors Negative jacobian etc Q 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 eg 4 then c Output element quantities e g stresses elseif isw eg 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
71. tructure of the routine UMESH1 is subroutine umesh1 uprt c User defined routine to input mesh data to FEAP implicit none include umaci h Contains UCT variable logical uprt Set name mesi to user defined if pcomp uct mesi 4 then uct xxxx Set user defined command name else c User execution function statements follow end if end 14 CHAPTER 4 USER FUNCTIONS 15 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 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
72. ture order for outputs Tla Mass interpolator a 0 Diagonal a 1 Consistent 8 q Loading intensity plates shells 9 To Stress free reference temperature 10 K 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 14ih 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 doo Material moduli 23 dag Material moduli 24 di Material moduli 25 dog Material moduli 26 d31 Material moduli 27 gi Material moduli 28 923 Material moduli 29 931 Material moduli 30 htype Heat flag Table 5 5 Material Parameters 39 CHAPTER 5 ADDING ELEMENTS Parameter Name Description 31 wv Orthotropic angle x principal axis 1 32 A Area cross section beam truss 33 T Inertia cross section beam truss 34 Iso Inertia cross section beam truss 35 ho Inertia cross section beam truss 80 J Polar inertia cross section beam truss 37 Ki Shear factor plate 38 Ko Shear factor plate 39 Non linear flag beam truss 40 Inelastic material model type 41 Yo I
73. ve with pAL pAL Mir Ma 3 Mp Mp1 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 residual 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 if isw eq 3 or isw eq 6 then c Compute element length L2 0 0d0 do i 1 ndm L2 12 x1 i 2 xl i 1 xx2 end do L sqrt L2 c Compute strain displacement matrix bb i 1 x1 1 2 xl i 1 eLr bb i 2 bb i 1 eps eps bb i 2 x 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 EA eps L body 0 5d0 L dm do i 1 ndm r i 1 body d 6 i bb i 1 sigA amp cmdxul 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 CH
74. y integer quantities and thus it was not deemed necessary to provide for integer history variables separately Although users may define new values for any 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 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 real 8 epl integer iepl neplts common ptdat6 ep1 20 iep1 2 20 neplts For problems in solid mechanics the linear momenta are stored as follows The linear momenta are comput
75. y 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 xl ndm nen 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 k j 1 Displacement ul ng j 2 Increment ult Un j 3 Increment u iy j 4 Rate ys j 5 Rate a J 6 Rate vn 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 per 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 29 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 select case isw case 1 c Input material parameters case default isw task Description Access Command 0 Output label SHOW ELE

FEAP - - A Finite Element Analysis Program

Contents

Download Pdf Manuals

Related Search

Related Contents