Home

FEAP Programmer Manual - Department of : Civil and

1. 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 Other values for use in writing material models can be obtained from parameters in common blocks For models which depend on position in the body the values of the reference and current coordinates for the constitutive point are passed in common block elcoor which contains the values in realx8 xref XCUr common elcoor xref 3 xcur 3 For models that may need to use an incremental formulation with Ae Epy En CHAPTER 4 USER FUNCTIONS 23 the array for strains may be dimensioned as eps 9 2 where the first 6 entries of eps 9 1 store the strains at t and the first 6 entries of eps 9 2 store those at tn The extra entries are not defined as they are provided only for use in the finite deformation form of the model described next Finite deformation models For finite deformation models the deformation gradient is passed to UMATLn in the ar gument array f 3 3 4 where f 3 3 1 defines Ru 3 3 2 defines F f 3 3 3 defines G 1 and f 3 3 4 defines G The deformation gradient is stored as f i J 1 Fig bata GJ 2 Etta f i J 3
2. Solution Array Names Numbers and Sized Element Array Names Numbers and Sizes Element connection array IX use for element e Element control array IE use for material number ma Element degree of freedom assignment array IEDOF use for material num perma x wk ee s ob b s Color Vable tor Plots eis Ste Se a es oe E Se es Arguments of FEAP Element Subprogram Task Options for FEAP Element Subprogram R Required O Optional E For eigensolutions FEAP common block definitions FEAP common block definitions a a a Material Parameters Material Parameters Material Parameters Material Parameters Material Parameters Tangent Parameters qu da an yad ab omba Say S Kub Element Plot Definition Subprograms Momenta and Energy Assignments Quadrature Tor triangles vs s usas oo s U S o Ge S il Quadrature for tetrahedra s Rob Ev Color pallet for FEAP plots bee a AR A R Values for control of plots S MGa oba Ka b k Ed Chapter 1 INTRODUCTION In this part of the FEAP manual some of the options to extend the capabilities of the program are d
3. 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 real 8 A ul ndf x1 ndm t1 s nst nst r nst if isw eq 0 and ior 1t 0 then c Output element description write Elmt 1 Element description Figure 5 1 FEAP Element Subprogram Part 1 32 CHAPTER 5 ADDING ELEMENTS 33 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 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 elseif isw eg 4 then c Output element quantities e g stresses elseif isw eq 5 then c Return Element mass matrix c s nst nst consistent matrix r ndf nen diagonal matrix elseif isw eq 6 then c Compute residual only c r ndf nen element residual elseif isw eq 7 then c Return Surface loading for element c s nst nst coefficient matrix c r ndf nst nodal forces elseif isw eq 8 then c Compute stress projections to nodes diagonal
4. 1 2 Uses of Common and Include Statements FEAP contains many COMMON statements that are used to pass parameters and small array values between subprograms For example access to the debugging parameter debug is facilitated through common debugs Users may either place the common statement as well as data typing statements directly in the routine or may use an include statement For debugging the statement would be 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 For example the common file name comblk h is defined as realx8 hr integer mr common comblk hr 1024 mr 1024 The arrays hr 1024 and mr 1024 serve to pass all dynamically allocated arrays be tween subprograms using a pointer array contained in the common array named np x lor for user defined arrays in up located in the include file pointer bh See Sec tion 3 for more details on use of pointers All include files are located in the directories values 1024 are necessary to ensure loops on arrays using pointers directly are considered as long CHAPTER 1 INTRODUCTION 4 include and include integer4 or include integer8 The subdirectory integer4 defines 32 bit integer parameters and integer8 defines 64 bit integer parameters It is highly recommended that users use include files ra
5. This prevents any unexpected execution of a statement that appears after the GO TO Some of the options for additional data passed through common blocks is shown in Figure 5 2 with each variable defined in Table 5 3 Also in Figure 5 3 the reference to common blocks using include statements is shown In the prototype routine the number of nodes on an element nen which is used to dimension ul is passed in the labeled common cdata Additional discussion is given below on use of some of the other data passed through the common blocks 5 1 Material property storage The material parameters to be stored in the array D with pointer np 25 may be input using the subprogram INMATE This subroutine is accessed by the statement CALL INMATE D TDOF NEV TYPE where D is the array storing the material parameters TDOF is returned as the parameter to access temperature NEV is the number of element history variables to allocate to NH1 and TYPE is the element type This routine inputs the commands as described in the user manual and stores the data for each material set into the D array elements as described in the following tables CHAPTER 5 ADDING ELEMENTS 37 Variable Definition o Page eject option head Title record numnp Number of mesh nodes numel Number of mesh elements nummat Number of material sets nen Maximum nodes element neq Number active equations ipr Real variable precision nstep Total number
6. UMACRn subprogram 29 Table of colors 31 Transient solution CTAN definition 47 Newton method 43 Non linear 37 Truss element Newmark method 51 Non linear theory 61 Tangent modulus 62 Theory 47 Variational equation 49 Weak form 49
7. available commands when the HELP command is used while in mesh solution or plot mode FEAP contains a large number of commands which are not commonly used by many users To control the default number of commands displayed to users the com mands have been separated into four levels 0 Basic 1 Intermediate 2 Advanced and 3 Expert The level to be displayed when using the HELP command is given may be set in the integer variable HLPLEV That is setting hlplev 1 Intermediate results in commands up to the intermediate level being displayed It is possible to raise or lower the level during execution using the command MANUal level where level is the numerical value desired When developing program modules it is often desirable to have output of specific quantities available e g tracking the change in some parameters during successive iterations FEAP provides for a switch to make the outputs active or inactive during an execution The switch is named debug and placed in integer ndebug logical debug common debugs ndebug debug CHAPTER 1 INTRODUCTION 3 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
8. i iste 4 Returns number projections end Figure 5 8 Element variable projection routine Chapter 6 UTILITY ROUTINES The FEAP system includes several subprograms that can assist developers in writing new modules In the next sections we describe some of the routines which perform numerical 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 6 2 3 Here only the description of subroutine calls is included together with the available options on number of points 6 1 1 One dimensional quadrature Line integrals may be evaluated using Gaussian quadrature in which the approximation to an integral is given as 1 L FE dE X gt KEW 6 1 1 1 1 where amp are quadrature points and W are the weights to be applied at each point The weights satisfy the condition S m 2 6 2 CHAPTER 6 UTILITY ROUTINES The Gauss Legendre formula has points which are all less than unity The subpro gram call CALL INTID L SW in which L is assigned an integer value between 1 and 5 returns the points and weights are returned in the two dimensional array SW 2 of type REAL 8 Points in SW 1 and weights in SW 2 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
9. quadratic and cubic order may be obtained using the subprogram call CALL TRISHP SS XL NDM IORD XJAC SHP where 3 3 6 5 1 2 1 4 2 3 Node Simplex 6 Node Element Figure 6 2 Triangular surface type elements in FEAP library CHAPTER 6 UTILITY ROUTINES 11 4 Node Element 3 4 10 9 3 11 8 8 6 12 7 1 5 2 1 5 6 2 8 Node Element 12 Node Element 4 7 3 8 6 1 5 2 9 Node Element 16 Node Element Figure 6 3 Quadrilateral surface type elements in FEAP library Parameter Description SS 3 Area coordinates Li Lo Ls XL NDM Element coordinates in local order NDM Spatial dimension mesh 2 or 3 IORD Order of interpolation 1 3 node 2 6 node 3 7 node 4 6 node 3 bubble 10 10 node cubic XJAC Jacobian transformation from x y to 7 SHP 3 NEL Shape functions and derivatives CHAPTER 6 UTILITY ROUTINES 78 The array SHP stores the values in the order SHP 1 i derivative with respect to or x SHP 2 i derivative with respect to 7 or y SHP 3 i shape function The parameter IORD defines the order of interpolation If it is 1 simple 3 node triangles with linear interpolation is 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 nodes 6 2 3 Shape functions in three dimensions The shape funtions for three dimensional tetrah
10. 6 UTILITY ROUTINES 82 Failure to include a pstyp statement will usually result in unpredictable plots of the mesh and contour values The known types of plots are 1 Point element with one node obtained by call CALL PLIPTI IEL 2 Line element with two nodes obtained by call CALL PLTLN2 IEL and for three node elements CALL PLTLN3 IEL 3 Triangular element with 3 nodes obtained by call CALL PLTRI3 IEL and for 6 nodes obtained by call CALL PLTRI6 IEL 4 Quadrilateral element with 4 nodes obtained by call CALL PLQUD4 IEL for 8 or 9 node elements the plot call is CALL PLQUD8 IEL and for 12 or 16 node quadrilaterals the call is CALL PLTQ16 IEL 5 Tetrahedral element with 4 nodes obtained by call CALL PLTET4 IEL CHAPTER 6 UTILITY ROUTINES 83 and for 10 node tetrahedra the call is CALL PLTET10 IEL 6 Brick element with 8 nodes obtained by call CALL PLBRK8 IEL and for 20 or 27 node bricks the call is CALL PLBRK27 IEL Using these and internal extraction of element surfaces the program is able to make some hidden surface plots in three dimensions 6 4 2 Element data plots Users may construct plots within their elements i e an ELMTnn and access using the plot command PLOT PELE vi v2 v3 In interactive mode in the plot environment it is only necessary to enter PELE vi v2 v3 The values entered in vi v2 v3 are optional and are passed to the element
11. Element coordinates in local order NDM Spatial dimension mesh 3 NEL Number of nodes on element 4 10 11 14 15 XJAC Jacobian transformation from xyz to Eng SHP 4 4 Shape functions and derivatives The array SHP stores the values in the same order as for the brick element 6 3 Eigenvalues for 3 x 3 matrix Three dimensional problems often require the solution of a 3 x 3 eigenproblem to generate principal values and directions FEAP includes a special routine to calculate the values and vectors for symmetric arrays The routine is used by a call to the subprogram as CALL EIG3 V D ROT On call to the routine V 3 3 is a REAL 8 array containing the symmetric array to be diagonalized On return the eigenvalues are contained in D 3 and the vectors for each value in the columns of the V array A Jacobi method is used with ROT an integer parameter returning the number of rotations to diagonalize The routine is quite efficient compared to any attempt to compute vectors after closed form solution of the cubic for roots In addition to the general eigensolution above FEAP includes options to compute principal values of a symmetric second order tensor for two and three dimensional problems In two dimensional use the call to CALL PSTR2D SIG PV is used where SIG 4 stores stresses in the order 011 022 033 012 and returns principal values and directions in PV 3 in the order or 02 and 6 where the angle is in degrees between
12. Table 6 1 Quadrature for triangles CALL TINT2D L LINT SW where L is a type indicator LINT returns the number of points and SW 4 is an array which returns three area coordinates and the quadrature weight SW 1 1 returns the area coordinate Lu as defined in 2 31 SW 2 1 returns the area coordinate La SW 3 1 returns the area coordinate La SW 4 1 returns the weight Table 6 1 describes the admissible types number and location of quadrature points 6 1 3 Three dimensional quadrature Three dimensional quadrature on brick domains may be performed by repeated one dimensional integration The three dimensional integrations are approximated by 1 L HE mQdednd 5 FE 6 4 m 1 1 where L is the total of all quadrature points A routine to compute n x n x n order Gauss Legendre quadrature is obtained by the call CALL INT3D L LINT SW where L is assigned to the number of points in each direction LINT is returned as the total number of points and SW 4 is an array containing the points and weights according to SW 1 1 contains values of the points SW 2 1 contains values of the CHAPTER 6 UTILITY ROUTINES 73 Type Number Location Points 1 1 Centroid O h 1 4 h 2 4 Interior O h 3 5 Interior O h Negative wt 4 11 Interior O h Negative wt 4 11 Nodal O h3 5 14 Interior O h 6 16 Interior O h 8 8
13. Truss The governing equations for a typical truss member element shown in Figure 5 4 are the balance of momentum equation Ao Os the strain displacement equation for small deformations Ab pAii 88 Os and a constitutive equation For example considering a linear elastic material the constitutive equation may be written as CHAPTER 5 ADDING ELEMENTS 48 Oss Ee 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 s is the coordinate along the truss member axis 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 u is a displacement in direction s e is an acceleration in direction s v css is a strain along the truss member axis and oss is the stress on a truss cross section The equations may also be deduced from the variational equation d d ou S t zeAnds Inn l L YL YL where Al contains the boundary and loading terms not associated with an element Where in addition to previously defined quantities we define e dis the spatial dimension of the truss 1 2 or 3 x are the Cartesian coordinates in the 4 directions L is the length of the truss member du is a virtual disp
14. Type 1 81 Method Type 2 82 Truss Rod quadrature number 83 Axial loading value 84 Constitutive start indicator 85 Polar angle indicator 86 Polar angle coord 1 87 Polar angle coord_2 88 Polar angle coord_3 89 Constitution transient type 90 41 Plane stress recovery 91 432 Plane stress recovery 92 a3 Plane stress recovery Table 5 7 Material Parameters 41 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 m Reference vector parameter 99 n3 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 Tso 129 Reference absolute temperature 130 nseg Number of hardening segments 131 148 Segment data sets epYisoHkin 149 Total variables on frame section 150 Piezoelectric flag 151 159 Piezoelectric data 160 Initial stress flag 161 166 cij Initial stresses constant 167 Tension compression only indicator 170 C Fung pseudo elastic model modulus 171 41 Fung model energy parameter 172 ag Fung model energy parameter 173 a3 Fung model energy parameter 174 ag Fung model energy parameter 175 as Fu
15. aslan 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 Ra see SA R Ba f R b ds gafos Ilan echt 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 o Ra 1 1 Oss A Ax pAL 2 1 Ui 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 os may be written us Al 1 Az HA AG 1 as For the linear elastic material a stiffness matrix may be expressed as where ae kij SECH Az NG E CHAPTER 5 ADDING ELEMENTS 51 The residual may now be written using a stiffness and mass matrix as d E 1 Jl kij uj m mil Vie Ri Fa H 3 E 7 5 2 Ma M23 Uiz with pAL pAL Tit m 5 Mig M21 3 6 For non linear material behavior the residual must be computed using Equation 5 1 with the stress replaced by the value computed from the constitutive equation The integration method for nodal quantities is taken as Newmark s method described in Section 5 2 The re
16. ix i end do i This avoids the need to specify appropriate fixed boundary conditions for the unused values Instead if one wishes to set the active degrees of freedom at each individual node of an element it is necessary to dimension the array as ix ndf x In this form the first column corresponds to the global pattern described above and columns 2 to nen 1 are associated with the local element nodes 1 to nen The element degrees of freedom are then assigned to each node individually by assigning a 1 for an active degree of freedom or O for an inactive one Note when using this option Do not make changes to the first column of the ix array Example 3 node element with 3 dof node Consider a problem with three degrees of freedom and three nodes on each element It is desired to have degrees of freedom 1 and 3 active on node 2 and degree of freedom 2 active on nodes 1 and 3 This is accomplished by setting the ix array values as ix 1 2 1 ix 2 2 1 ix 3 2 ix 1 3 1 2 ix 2 3 ix 3 3 1 ix 1 4 I For node 3 ix 2 4 1 ix 3 4 0 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 Setting element plot sequence Utility routines are also provided to provide the necessary list of nodes needed to properly draw the mesh for each element type during plot outputs The nam
17. projection weight wt nen c s nen projection values st nen c 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 De 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 34 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 labeled as R shown above If elements have local variables that need to be retained between subsequent time steps history variables may be defined as described in Section 5 6 In this case it is necessary to code task 12 if special transformations of the variables are required otherwise merely return with no changes and if any of the parameters have non zero initial values task 14 is used to set these values zero values are set by default Finally if special plotting options are desired it may be necessary to program task 20 note that contours for element variables such as st
18. started over On the other hand if convergence occurs during the time step and the value of rmeas is smaller than 1 25 the time step is adjusted according to A 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 Atreo AES 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 CHAPTER 4 USER FUNCTIONS 28 subroutine umanil c User defined routine to manipulate mesh data for FEAP implicit none include umaci h 1 Contains UCT variable 6 Set name 1 to user defined if pcomp uct mani 4 then uct Set user defined command name else User execution function statements follow end if end Figure 4 6 Sample UMANLn module where n has a value between 0 and 9 must be written compiled and link
19. the statement rotyp mm in the section of ELMTnn where isw 1 This parameter is located in the common erotas which has the structure realx8 Xin realx8 rots rvel racc thkl CHAPTER 5 ADDING ELEMENTS 58 integer rotyp common erotas x1n 9 9 4 amp rots 3 9 2 rvel 3 9 2 racc 3 9 2 thk1 9 rotyp The other entries in the common are arrays that return values for each element to treat the rotation values and rates We shall return to their description after describing the treatment of the global nodal data for rotations 5 7 1 Nodal rotation treatment UROTmm subprogram The nodal rotation data is stored in the array xlg which is dimensioned as x1g 9 6 nump For node ng the entries in x1g are stored as follows Component I O Description XLG 1 9 1 ng I Rotation matrix A at time ty Alterna tively entries 1 to 4 may be used to store a quaternion Rotation matrix A at time tn Rotation matrix at time tn 1 Rotation increment angle A Rotation rate w at time tn Rotation acceleration a at time tn Rotation angle O Rotation rate d n at time tara Rotation acceleration Gra at time tara Rotation matrix A at time to XLG 1 9 2 ng XLG 1 9 3 ng XLG 1 3 4 ng XLG 4 6 4 ng XLG 7 9 4 ng XLG 1 3 5 ng XLG 4 6 5 ng XLG 7 9 5 ng XLG 1 9 6 ng While storage is provided for the 3 x 3 rotation matrices the representation may also be specified
20. 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 CHAPTER 6 UTILITY ROUTINES 84 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 fil 2 Move to point 3 Draw to point Table 6 4 Values for control of plots 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 21 1 21 and 22 y 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 a panel is drawn with stand
21. to zero Thus if the array is to be used only once it need not be set to zero before accumulating additional values If the array is to be reused or resized see below it must be reinitialized prior to accumulating any additional values Use of these subprograms controls the assignment of memory space for all arrays such that no CHAPTER 8 ALLOCATING ARRAYS 10 conflicts occur between hr and mr referenced arrays Each routine which makes direct reference to an allocated array using a pointer e g hr np 43 or mr up 1 must contain include files as include pointer h include comblk h As an example for the use of the above allocation scheme consider a case where it is desired to allocate a real double precision array with length NUMNP number of nodes in mesh and an integer array with length NUMEL number of elements in mesh The parameters NUMNP and NUMEL are contained in COMMON CDATA and available using the include file cdata h The new arrays are defined using the temporary names TEMP1 and TEMP2 which have numerical locations 111 and 112 respectively The two arrays are allocated using the statements setvar setvar palloc 111 1 nump 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 th
22. x and the 1 axis This routine does not use SIG 3 CHAPTER 6 UTILITY ROUTINES 81 In three dimensions the principal values are obtained using the call CALL PSTR3D SIG PV where SIG 6 stores stresses in the order 011 022 033 012 023 031 and returns principal values in PV 3 in the order 01 02 03 Roots are ordered from most positive to most negative 6 4 Plot routines Several options exist in the FEAP system to create graphical plots for data and results 6 4 1 Mesh plots FEAP has plot capabilities to represent some standard element shapes For continuum elements where the shape of the element is identical to the space dimension of the mesh i e NDM it is not necessary to provide any extra descriptions However if the dimension of the element topology is different from the mesh dimension it is necessary to add the include statement include eldata h and the statement pstyp pdim within the ISW 1 part of the element routine Here pdim is the dimension of the element topology For example in a three dimensional shell problem where NDM 3 and the element topology is two dimensional the statement is given as pstyp 2 Provided the nodal numbering of an element is as described in the FEAP User manual i e numbered with vertex nodes first followed by mid side nodes then face nodes and finally internal nodes the program can use the actual number of nodes on the element to draw each element CHAPTER
23. 1 9 nl 1 xlg 1 9 1 ng xln 1 9 n1 2 xlg 1 9 2 ng xln 1 9 n1 3 xlg 1 9 3 ng xln 1 9 n1 4 xlg 1 9 6 ng rots 1 3 nl 1 xlg 1 3 4 ng rots 1 3 nl 2 xlg 1 3 5 ng rvel 1 3 nl 1 xlg 4 6 4 ng rvel 1 3 nl 2 xlg 4 6 5 ng racc 1 3 nl 1 xlg 7 9 4 ng racc 1 3 n1 2 xlg 7 9 5 ng where n1 is a local node number between 1 and 9 the maximum provided in the current erotas and ng is the global node number associated with each local number CHAPTER 5 ADDING ELEMENTS 60 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 12 Momenta and Energy Assignments Using the above data structure one can program the updates in any manner that does not conflict with the time treatment The only interface to FEAP is through how the increment du 4 6 n is defined 5 8 Energy Computation FEAP elements provide an option to accumulate the total momenta and energy during the solution process The values are accumulated in the array EPL 20 when the switch parameter isw is 13 and written to a file named Pxxxx ene where xxxx is extracted from the problem input filename whenever the solution command TIME is used The array EPL 2 is in the common block named ptdat6 which has the structure real 8 epl integer iepl neplts common ptdat6 ep1 20 0 iep1 2 200 neplts Fo
24. 1 with the remainder in the interior of the interval A routine to perform quadrature is obtain by using the call CALL INTIDL L SW in which L is assigned an integer value between 1 and 6 The values of the points and weights are returned in the two dimensional array SW Points in SW 1 and weights in SW 2 6 1 2 Two dimensional quadrature Two dimensional quadrature on quadrilateral domains may be performed by repeated one dimensional integration The two dimensional integrations are approximated by ff temagan 5 Kem W 6 3 where L is the total of all quadrature points A routine to compute n x n order Gauss Legendre quadrature is obtained by the call CALL INT2D L LINT SW where L is assigned to the number of points in each direction LINT is returned as the total number of points and SW 3 is an array containing the points and weights according to SW 1 1 contains values of the points SW 2 1 contains values of the points 7 and SW 3 1 contains values of the weights W Two dimensional quadrature on triangles may be performed using the subprograms call CHAPTER 6 UTILITY ROUTINES 72 Type Number Location Points 1 1 Centroid O h 3 3 Mid sides O h 3 3 Interior O h 4 4 Interior O h Negative Wt 6 6 Nodal O h3 6 6 Interior O h F 7 Interior O h 7 7 Nodal O h9 12 12 Interior O h 13 13 Interior O h3 Negative Wt
25. ADDING ELEMENTS 57 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 De solution command TIME These access the task options isw equal to 3 or 6 and 14 respectively If a user adds a new option for which it is desired to save the history variables it is necessary to set the variables hflgu and h3flgu to true as required if no update is wanted the variables should be set to false These parameters are located in logical hflgu h3flgu common hdatam hflgu h3flgu 5 7 Elements with Finite Rotation Parameters When considering structural elements that undergo large displacements it is usually necessary to treat the rotation parameters for large angle changes The nodal pa rameters for this case are a rotation vector 6 and the finite rotations are given as an orthogonal matrix A A Any exp 6 A in which 6 denotes a skew matrix given as o o 83 o 0 0 0 9 6 9 0 The actual method used to update the rotations and their increments must be specified when writing the element module ELMTnn and is performed by a user subprogram named UROTmm where mm is a number between 01 and 10 To specify which routine is to be used it is necessary to include
26. CON use 18 UMATIn subprogram 20 UMATLn subprogram 20 Material property variables 36 Mesh array names 11 Mesh input UMESHn subprogram 15 Mesh manipulation UMANIn subprogram 27 Mesh plots 81 Brick element 83 Line element 82 Point element 82 Quadrilateral element 82 Tetrahedral element 82 Triangular element 82 Newton solution 90 INDEX Element residual array storage 45 Element tangent array storage 45 Non linear problems 43 Tangent array definition 44 Numerical integration 70 One dimension 70 Three dimension 72 Two dimension 71 Partial list of element common block vari ables 64 Plot command UPLOT subprogram 85 UPLOTn subprogram 29 Plot utility function Fills 84 Lines 84 Plot utility functions Color 30 Fills 29 Lines 29 Precision Integer 1 Real 1 Principal stress Three dimensional 80 Two dimensional 80 Problem size 1 Quadrature 70 One dimension 70 Three dimension 72 Two dimension 71 Setting colors 30 Setting options Help level 2 Parsing statements 1 Plot prompts 2 PostScript 2 91 Shape functions 73 Brick isoparametric 78 One dimensional cubic Hermitian 75 One dimensional isoparametric 74 One dimensional natural derivative 74 Quadrilateral isoparametric 75 Tetrahedral isoparametric 79 Three dimensional isoparametric 78 Triangular isoparametric 76 Two dimensional isoparametric 75 Solution array names 11 Solution command
27. EMENTS 64 character 4 o head common bdata o head 20 integer numnp numel nummat nen neq ipr common cdata numnp numel nummat nen neq ipr integer nstep niter nform naugm titer taugm tform common counts nstep niter nform naugm titer taugm tform integer iaugm iform intvc iautl nstepa nsplt common counts iaugm iform intvc iautl nstepa nsplt realx8 dm integer n ma mct iel nel pstyp eltyp common eldata dm n ma mct iel nel pstyp eltyp real 8 tt common elplot tt 1000 real 8 bpr ctan psil common eltran bpr 3 ctan 3 psil real 8 ut common eluser ut 1000 integer nhi1 nh2 nh3 hti1 ht2 ht3 1 int 4 or int 8 common hdata nh1 nh2 nh3 ht1 ht2 ht3 int 4 or int 8 integer nlm plm nge common hdata nlm plm nge integer ior iov ilg common iofile ior iov ilg logical Keepfl common iofile keepfl integer nph ner int 4 or int 8 real 8 erav int jshft common prstrs nph ner erav j int 3 jshft integer ndf ndm nen1 nst nneq ndl nnlm nadd common sdata ndf ndm nen1 nst nneq ndl nnlm nadd realx8 ttim dt c1 c2 c3 c4 cb chi dtcr integer idyn0 common tdata ttim dt c1 c2 c3 c4 c5 chi dtcr idyn0 integer n Up common pointer np 400 up 200 real 8 hr integer mr common comblk hr 1024 mr 1024 Figure 5 2 FEAP Element Common Blocks CHAPTER 5 ADDING ELEMENTS include bdata h include cdata h include counts h include eldata h include
28. FEAP A Finite Element Analysis Program Version 8 4 Programmer Manual Robert L Taylor Department of Civil and Environmental Engineering University of California at Berkeley Berkeley California 94720 1710 E Mail rltQce berkeley edu Revised June 2014 Contents Introduction 1 1 Setting Program Options ebe ARA E a 1 2 Uses of Common and Include Statements Data Input and Output 2 1 Parameters and Expressions 2 sd dary a a eae GG 22 Array Outp ts see rana an s s d n Allocating Arrays User Functions 4 1 Mesh Input Functions UMESHn 4 1 1 Command line TX data 4 2 User Material Models 4 2 1 The UMATIn Module 4 2 2 The UMATLn Module 0000000004 4 2 3 Auto time step control a as fon NALANG EE 4 3 Mesh Manipulation Functions UMANIn 4 4 Solution Command Functions UMACRn 4 5 Plot Command Functions UPL Tn Adding Elements 5 1 Material property storage suco a gt Epa a ER ER pd 5 2 Non linear Transient Solution Forms 5 3 Example 2 Node Truss Element AREA A 5 3 1 Theory for a Truss EECHER 5 4 Additional Options in Elements oa Task Options lt eoe s a BRAND A Ha Task O ETH qu e E 7 7 5 5 Projection of element variables to nodes 5 6 Elements with His
29. Giz tnga dy and NAIJA Gist Killa is where G y are displacement gradients Stress and moduli are to be returned in the ar gument arrays dimensioned as sig 6 and dd 6 6 respectively Stresses and moduli are returned using Voigt notation where stresses are ordered as T o o 022 033 912 023 031 with corresponding order for the moduli All real values are in 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 deformation gradient are as mentioned above passed in the array f 3 3 4 and the determinant of the deformation gradient in the parameter theta 4 where 04 detF y1 and 0 det Fn In addition 0 det Fn 41 1 and 04 detF 1 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 Other values for use in writing material models can be obtained from parameters in common blocks For models which depend on position in the body the values of the CHAPTER 4 USER FUNCTIONS 24 subroutine umatl1 eps theta td d ud hn hi nh ii istrt sig dd isw Cranes c Purpose User Const
30. HAPTER 2 DATA INPUT AND OUTPUT 6 input file unit number if negative input from keyboard iow output file unit number ilg solution log file unit number If an error occurs during input from the keyboard FEAP returns a value of true for the function and a user may reinput the record if the implied loop shown above is used For inputs from a file the program will stop and an error message indicating the type of error occurring and the location in an input file is written to the output file The input routines return data in a real 8 array td If any td i is to be used as an integer or real 4 quantity it must be cast to the correct type That is the following operations should be used to properly cast the variable type 1 4 t real 8 td 5 integer j logical errck pinput errck pinput td 5 nint td 1 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 1 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 rea
31. MESHn To add a mesh input command a subprogram with the name UMESHn where n has a value between 0 and 9 must be written compiled and linked with the program The basic structure of the routine UMESH1 is The parameter TX is a character array which is assigned by the input and 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 to assign the name of the command The default name is 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 1 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 15 CHAPTER 4 USER FUNCTIONS 16 subroutine umeshi tx uprt c c Purpose User defined routine to input mesh data to FEAP c Inputs c tx Co
32. Node amp Face O h Table 6 2 Quadrature for tetrahedra points 7 and SW 3 1 contains values of the points and SW 4 1 contains values of the weights W Three dimensional quadrature on tetrahedra may be performed using the subprograms call CALL TINT3D L LINT SW where L is a type indicator LINT returns the number of points and SV 5 4 is an array which returns three area coordinates and the quadrature weight SW 1 1 returns the volume coordinate L as defined in 2 3 SW 2 1 returns the volume coordinate Lou SW 3 1 returns the volume coordinate La SW 4 1 returns the volume coordinate Lai SW 5 1 returns the weight W Table 6 2 describes the admissible types number and location of quadrature points 6 2 Shape function subprograms Finite element approximations commonly use shape function subprograms to perform computations of the functions and their derivatives at preselected points often the quadrature points FEAP includes options to obtain the shape functions for some low order elements linear and quadratic order in one and two dimensions and linear shape functions for three dimensions In addition a cubic Hermitian interpolation routine is available The calling arguments for routines is summarized below CHAPTER 6 UTILITY ROUTINES 74 6 2 1 Shape functions in one dimension The shape funtions for one dimensional elements as shown in Fig 6 1 may be com puted using the shape function r
33. SH1 the first argument GETData must match the name assigned to UCT and will also be in TX 1 The second parameter will be in TX 2 and the third in TX 3 To recover the numerical value for the third parameter the statement statements real 8 ctl call setval tx 3 15 ctl may be used to assign the real value 35 0d0 to ctl If necessary the real value for 61 can be cast into an integer using itl nint ctl If more than 8 items are desired on the input line it is possible to recover their values from the character string yyy 256 which has been parsed into columns with width 15 characters 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 To recover their values the statements CHAPTER 4 USER FUNCTIONS 18 include chdata h character lct 15 x15 real 8 rt1 15 integer itl are added to the user function and the items recovered in the else option of the function using the statements 1ct 1 yyy 16 30 call setval yyy 31 45 15 rt1 1 would assign 1ct 1 values from the second set of 15 characters and rt1 1 to the third set of 15 characters In this case 1ct 1 tx 2 and rtl 1 would have the same value as ctl above If users wish to add more than 10 material models it is possible to use the user function UMESH which has the form 4 2 User Material Models Users may add material models to elements by appending subprograms UMA
34. SHOW ELEM 1 R Input d parameters Mesh MATE n 2 O Check elements Soln CHECk 3 R Compute tangent residual Soln TANG Store in S r UTAN 4 O Output element variables Soln STRE 5 E Compute cons lump mass Soln MASS Store in S r MASS LUMP 6 R Compute residual Soln FORM REAC Plot REAC 7 O Surface load tangents Mesh SLOAd 8 O Nodal projections Soln STRE NODE Plot STRE PSTR 9 O Damping Soln DAMP 10 O Augmented Lagrangian update Soln AUGM 11 O Error estimator Soln ERRO 12 R History update Soln TIME For special history treatments else return 13 O Energy momentum Soln TPLO ENER 14 R Initialize history BATCh INTEr 15 O Body force Mesh BODY 16 O Jintegrals Soln JINT 17 O Set after activation Soln ACTI 18 O Set after deactivation Soln DEAC 19 NOT AVAILABLE used in modal base BASE Uses isw 5 in element 20 O Element plotting Plot PELM 23 O Compute element loads only ARCL 25 O Zienkiewicz Zhu projection Soln ZZHU 26 R Used to compute mesh boundary Called by default Table 5 2 Task Options for FEAP Element Subprogram R Required O Optional E For eigensolutions CHAPTER 5 ADDING ELEMENTS 36 select case isw case 1 c Input material parameters case default end select the CASE DEFAULT should not perform any operations unless all options are programmed Finally if the form go to 1 2 isw return is used the RETURN statement should always be included as shown
35. TIn 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 viv2 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 CHAPTER 4 USER FUNCTIONS 19 C C C C o o logical function umesh cc tx prt Ene a me oth 03 9 i Purpose User mesh command interface Inputs CC User command option tx x Command line input data prt Output if true Outputs none Data stored by user development Saan Kees implicit none logical prt pcomp character 4 tx 15 Match on USER Add as many checks as desired with user if pcomp cc xxxx 4 then Provide name for xxxx umesh true Activate command elsei
36. ard line drawing commands starting with CALL PLOTL X1 Y1 Z1 1 CHAPTER 6 UTILITY ROUTINES 89 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 i Yi zi specification is made In this case a user may wish to pass the value of some solution variable as the z value scaled so that it will make sense relative to the x 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 parameter i e variable CS is available in labeled common PVIEW In this case one can add the statements assuming here that the displacements correspond to the coordinate directions DO NE 1 NEL DO I 1 NDM XP I NE XL I NE CS UL I NE END DO I END DO NE NEL is the number of connected nodes to each element and is passed through labeled common ELDATA before performing any deformed plots and then plot the appropriate values of XP Indeed this may always be performed as the valu
37. cify new name for model c Input output user data and save in ud array else c Set values of nl if required ni ml write iow User Constitutive Inputs E vv 1 ud 1 vv 1 Parameter from input on command name endif end Figure 4 3 Sample UMATI1 module CHAPTER 4 USER FUNCTIONS 22 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 Small deformation models For small deformation models the strains are passed to UMATLn in the argument array eps 6 and stored in the order T e c 22 33 712 V23 oa where yj 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 real values are in 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 the trace of the strain in the parameter theta Thus 0 64 11
38. df 11 do ii 1 2 11 0 1 2 do i 1 ndm do 1 ndm s itil j j1 db i ii bb j jj end do end do ji ji ndf end do 11 il ndf end do c Compute mass terms and correct for inertial effects cmd ctan 3 rhoA L 3 0d0 cmo cmd 0 5d0 do i 1 ndm j i ndf 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 6 end do endif Figure 5 7 Truss Tangent Matrix Part 2 CHAPTER 5 ADDING ELEMENTS 69 subroutine slcn2d sig shp xsj sg lint nel nes p s gr E nunun sian nuna es mena 1 c Purpose Project element variables to nodes c Inputs sig nes Stresses at quadrature points c shp nel Shape functions at quadrature points xsj x Volume element at quadrature points sg 3 Gauss points 1 2 and weights 3 c lint Number of quadrature points c nel Number nodes on element c nes Dimension of stress array c Outputs c p nen Weights for lumped projection c s nen Integral of variables Go EE implicit none include cdata h 1 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 pi xg s i 1 s i 1 sig 1 1 xg 8 1 2 s i 2 sig 2 1 xe 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
39. e whereas if the array is not found it is false The precision will be either one 1 or two 2 for INTEGER or double precision REAL 8 quantities respectively Thus the above coordinate call will return xpre as 2 and xlen will be the product of the space dimension of the mesh and the total number of nodes in the mesh The first coordinate z may be given as x1 hr xpoint any other coordinates at nodes may also be recovered by a correct positioning in later words of hr For example y is located at hr xpoint 1 The use of pgetd can lead 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 U
40. e 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 CHAPTER 6 UTILITY ROUTINES 86 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 Chapter 7 Adding a user solver ADDING USER SOLVERS There are several public domain linear equation solution routines available at various internet locations Examples are SuperLU umfpack xxx to name three To access any of these solvers it is necessary to add user modules named umacr1 f and usolve f to FEAP The module umacrx f x ranges between 0 and 9 has the basic form subroutine umacr1 l1ct ctl prt include setups h for parameter solver include umaci h for parameter uct logical prt character lct 15 real 8 ct1 3 if pcomp uct mac1 4 then uct Set name of command for solver else if pcom
41. e 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 the 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 UMATIn Module A sample for the UMATL1 module with arguments defined for small deformation is shown in Fig 4 4 and for arguments defined for finite deformation in Fig 4 5 This subpro gram will be called by many of the elements included within FEAP if a user model CHAPTER 4 USER FUNCTIONS 21 subroutine umatil type vv d ud n1 n3 g Bee io SS apos ooo 1 c Purpose User material model interface c Inputs c type Name of constitutive model character variable c vv Parameters user parameters from command line c d Program material parameter data c Outputs c ni Number history terms nhi nh2 c n3 Number history terms nh3 c ud User material parameters o eebe implicit none include iofile h logical pcomp character type 15 integer ni n3 real 8 vv 5 d ud c Specify type of user model if pcomp type mat1 4 then type E 1d Spe
42. e with the name UMACRn where n ranges from 0 to 9 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 The parameters CTL 3 are used to pass the three parameter values read respectively Again the name xxxx should be selected to not conflict with existing plot command names and will appear whenever HELP is issued CHAPTER 4 USER FUNCTIONS 30 subroutine uplot0 ct1 c User plot command function implicit none include umaci h Contains the variable UCT real 8 ct1 3 c Set command word if bcomp uct p1t0 4 then uct else c User plot command statements are placed here endif end Figure 4 8 Sample UPLOTn module 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 1 to 3 1 starts a filled panel 2 dra
43. ed with the program The basic structure of the routine UMANT1 is 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 1 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 CHAPTER 4 USER FUNCTIONS 29 subroutine umacr0 1ct ct1 c User solution command function implicit none include umac1 hb Contains the variable UCT character lct 15 real 8 ct1 3 Set command word if pcomp uct macO 4 then uct else c User command statements are placed here endif end Figure 4 7 Sample UMACRn module 4 4 Solution Command Functions UMACRn In a similar manner users may add solution commands to the program by adding a routin
44. edral elements as shown in Fig 6 4 and brick elements as shown in Fig 6 5 may be computed using the shape function routines described below 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 2 4 Node Simplex 10 Node Element Figure 6 4 Tetrahedron solid type elements in FEAP library CHAPTER 6 UTILITY ROUTINES 79 20 Node Element 27 Node Element Figure 6 5 Brick solid type elements in FEAP library 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 En The array SHP stores the values in the order SHP 1 i derivative with respect to 2 SHP 2 i derivative with respect to y SHP 3 i derivative with respect to z SHP 4 i shape function Three dimensional Co isoparametric interpolation on tetrahedra of linear order i e 4 node elements may be obtained using the subprogram call CHAPTER 6 UTILITY ROUTINES 80 CALL TETSHP SS XL NDM NEL XJAC SHP where Parameter Description ss 4 Volume coordinates Ly Ls L3 Ly XL NDM 4
45. eir pointers e g in a subroutine call as include pointer h include comblk h call subname hr np 111 mr np 112 Note the use of hr and mr for the double precision and integer references respectively Also the use of the pointers avoids a need to include the array reference until it is needed in a computation A short list of the mesh arrays available in FEAP is given in Table 3 1 for solution arrays in Table 3 2 and for element arrays in Table 3 3 The names of all active arrays in any analysis may be obtained using the SHOW DICTionary solution command The array IX nen1 numel is used to store basic information for each element in the mesh related to the nodal connections and material data requirements In addition See the subprogram palloc f in the program directory for the names and numbers of existing arrays CHAPTER 8 ALLOCATING ARRAYS 11 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 1 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 di
46. elplot h include eltran h include hdata h include iofile h include prstrs h include tdata h include pointer h include comblk h Figure 5 3 FEAP Element Common Blocks using Includes Figure 5 4 2 Node Truss Element CHAPTER 5 ADDING ELEMENTS 66 if isw eq 3 or isw eq 6 then c Compute element length L2 0 0d0 do i 1 ndm L2 L2 1 4 2 x1 i 1 2 end do L sqrt L2 c Compute strain displacement matrix bb i 1 x1 i 2 xl1 i 1 Lr bb i 2 bb i 1 eps eps bb i 2 ul i 2 1 ul i 1 1 end do c Compute mass terms cmd rhoA L 3 0d0 cmo cmd 0 5d0 Form body inertia force vector dm prop ld sigA EAxepsxL body 0 5d0 L dm do i 1 ndm r i 1 bodyxd 6 i bb i 1 sigA amp cmdxul i 1 5 cmoxul i 2 5 r i 2 bodyxd 6 i bb i 2 sigA amp cmoxul i 1 5 cmd ul i 2 5 end do Figure 5 5 Element residual for two node truss CHAPTER 5 ADDING ELEMENTS if isw eq 3 then c Compute element length L2 0 0d0 do i 1 ndm L2 L2 xl i 2 x1 i 1 xx2 end do L sqrt 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 1 4 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 67 CHAPTER 5 ADDING ELEMENTS 68 Compute stiffness terms ndm 4 or n
47. en 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 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 d OR i OR Ov OR d 7 Ov 7 Note that from the Newmark formulas 1 V Ov Y du bAt du da du PB At in which z is the Kronnecker delta identity matrix for the k j nodal pair From the residual we observe that OR OR OR Ky e 8 5 Mg Ov J da define the tangent stiffness damping and mass respectively Thus for the Newmark algorithm the total tangent matrix in terms of the incremental displacements is 1 Cu 4 5 Ki 4 BAP Y BAt CHAPTER 5 ADDING ELEMENTS 45 For other choices of increments the tangent may be written in the general form Sij CG K C2 Cy C3 M where the c are scalar quantities involving the integration parameters of the method selected and At Thus any one step integrator may be considered and will affect only the specification of the constants in the tangent In FEAP the element tangent matrix S is stored as a two dimensional a
48. ept 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 8 hr integer mr common comblk hr 1024 mr 1024 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 4 GByte range The pointers in the integer8 subdirectory have 64 bit lengths and should be used for additional index ranges The arrays hr and mr are used to establish addresses only and not to physically store data This mechanism permits references to elements in arrays which have positions relative to hr or mr that may be after or before 1 Thus FEAP must be compiled without strict array bound checking Size of problems is limited only by the available memory in the computer used CHAPTER 8 ALLOCATING ARRAYS 9 When using 64 bit pointers users must be careful to always define the addre
49. es of the routines are listed in Table 5 11 and each routine is called as CHAPTER 5 ADDING ELEMENTS 53 Routine Name Description PLTLN2 2 node line element PLTLN3 3 node line element PLTRI3 3 node triangular element PLTRI6 6 node triangular element PLQUD4 4 node quadrilateral element PLQUD8 8 or 9 node quadrilateral element PLTQ16 16 node quadrilateral element PLTET4 4 node tetrahedron element PLTET10 4 node tetrahedron element PLBRK8 8 node brick element PLBK27 27 node brick element PLBKPQR 64 node brick element Table 5 11 Element Plot Definition Subprograms call plname iel where iel is an integer parameter defined in common eldata If no call to a sub program is included each element is assumed to be a 4 to 9 node quadrilateral and a 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 11 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 real 8 tt common elplot tt 1000 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 CHAPTER 5 ADDING ELEMENTS 54 tt 1 tt 2 EA eps Eleme
50. es 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 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 While this wastes storage for integer variables experience indicates there is little need to save many integer quantities and thus it was not deemed necessary to provide for integer history variables separately CHAPTER 5
51. escribed 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 feap84 f module controls the specification of the ratio of REAL to INTEGER variables For typical UNIX and PC systems all real variables should be twice as large as integers and IPR is set to 2 For systems in which INTEGER 8 variables are used set by compiler option the IPR parameter is set to 1 Any error in setting this parameter may lead to incorrect behavior of the program consequently do not reset the parameter unless a careful assessment of compiler behavior has been made 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 da
52. f ikasa endif end Figure 4 2 Sample UMESH module CHAPTER 4 USER FUNCTIONS 20 MATErial ma SOLId ELAStic IS Tropic e nu UCONstitutive xxxx vi v2 No standard commands should follow the UCON command Alternatively users may input elastic properties as part of their UMATIn module If the user routine does input additional data records after the UCON record and these are terminated by a blank record a second blank record will be needed to discontinue material data input for this set In all cases at least one blank record is always needed to terminate the input of standard options for the material set 4 2 1 The UMATIn Module A sample module for a user constitutive model is shown in Fig 4 3 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 fifth argument of UMATIn In the current version there are 150 words of double precision values available by default Additional values may be allocated by assigning a larger valu
53. hod The Basis vol ume 1 Butterworth Heinemann Oxford 5 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 6 M Gellert and R Harbord Moderate degree cubature formulas for 3 D tetrahedral finite element approximations Communications in Applied Numerical Methods 7 487 495 1991 89 Index Adding elements 32 ELMTnn subprogram 32 Arguments on subprogram ELMTnn 34 Array allocation 8 PALLOC subprogram 9 Removing arrays 10 UALLOC subprogram 9 Color in plots 30 Common block definitions 38 Common statements 3 Data input Converting text to real 7 PINPUT 5 TINPUT 6 Data output Integer arrays 7 Real arrays 7 Debug output 2 Eigenvalues for 3 x 3 matrix 80 Element array names 12 Element connection array 12 Element control array 13 Element degree of freedom array 13 Element energy computation 60 Element history variables 55 Accessing for each element 56 Element options Projection to nodes 54 Task 1 ISW 1 51 Task 6 ISW 6 53 TPLOt definitions 53 Element plots User 83 Element task definitions 35 Finite rotation 57 Include statements 3 Locating arrays PGETD 14 Material models 18 Auto time step control 25 History variables 25 Internal variables 25 U
54. in terms of quaternions for which only 4 components are necessary In this case the 9 entries may be divided into two 4 entry quantities if required Indeed the space may be used in anyway necessary provided no conflict in the way each time value is associated to the data Note that sufficient storage is available to define integration methods for which the rotation is defined at an intermediate time tn a For a typical node n in the mesh the location of the entries in the x1g array are obtained from CHAPTER 5 ADDING ELEMENTS 59 ng mropt n 2 and the routine UROTmm is called as call urotmm xlg 1 1 ng xlg 1 2 ng xlg 1 3 ng amp xlg 1 4 ng xlg 1 5 ng amp x1g 4 4 ng xlg 4 5 ng amp xlg 7 4 ng xlg 7 5 ng du tsw where du 1 3 are the solution increments for rotation from the solver and tsw is the time update switch which is set according to tsw 1 Initialize for new time step tsw 2 Update within a time step tsw 3 Back up to beginning of time step The entry u 1 is the location for the first entry in the rotation vector 0 5 7 2 Local nodal rotation treatment When each element that is associated with nodal rotation parameters is called the rotation data is transferred to local storage in a manner similar to treatment of trans lations The local data is passed to each element using the common erotas defined above The entries in the local arrays are extracted from the global array according to xIn
55. inear 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 R RY A CIR dall 0 4 4 da 3 where a is one of the variables at time t 1 e g Uj tn41 We define e _2R iw 4 da and solve SO do RO ij j H The solution is updated using RAI o a CHAPTER 5 ADDING ELEMENTS 44 In the above k is the iteration number for the Newton algorithm To start the solution for each step FEAP sets oi ia alta where a quantity without the k superscript represents a converged value For a linear problem Newton s method converges in one iteration Computing the residual after one iteration must yield a zero value to within the roundoff of the computer used For non linear problems a properly implemented Newton s method must exhibit a quadratic asymptotic rate of convergence Failure of the above performance for linear and non linear cases implies a programming error in an implementation or lack of a consistently linearized algorithm i e S is not an exact derivative of the residual In a non linear problem Newmark s method may be parameterized in terms of in crements of displacement velocity or acceleration From the Newmark formulas the relations du da and dv yAtda define the relationships between the increments Note that only scalar multipliers involving 0 y and At are involved betwe
56. 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 x1 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 to add these terms The programming of the tangent array however must distinguish between cases in which transient e g inertial loads are present and those in which they are omitted The different cases are implemented in FEAP by making appropri ate assignments to the c parameters To facilitate the programming of the tangent array 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 10 Thus in solid mechanics applications the tangent matrix is defined in an element
57. ion a 0 Diagonal a 1 Consistent 814 Loading intensity plates shells 91 Stress free reference temperature 10 lk Shear factor plates shells beams 11 by Body force volume in 1 directions 12 bo Body force volume in 2 directions 13 bs Body force volume in 3 directions 14 h 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 da Material moduli 24 di Material moduli 25 dos Material moduli 26 d3 Material moduli 27 gi Material moduli 28 993 Material moduli 29 gan Material moduli 30 htype Heat flag Table 5 5 Material Parameters 39 CHAPTER 5 ADDING ELEMENTS Parameter Name Description 31 4 Orthotropic angle x principal axis 1 32 A Area cross section beam truss 33 ly Inertia cross section beam truss 34 Inertia cross section beam truss 35 h Inertia cross section beam truss 307 al Polar inertia cross section beam truss 37 Ki Shear factor plate 38 Shear factor plate 39 Non linear flag beam truss 40 Inelastic material model type 41 Yo Initial yield stress Mises 42 Wns Final
58. ion operations for each element under ISW 8 These are performed locally for 3An implementation of the Zienkiewicz Zhu projection method is implemented using ISW 24 CHAPTER 5 ADDING ELEMENTS 59 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 which are not associated to any particular time All history variables are associated with the allocation name H which has a pointer value 49 Users are not permitted direct access to the data stored as H of course it is possible to access from hr np 49 but this should not normally be attempted Before calling the element routine for each element FEAP transfers the required history variable to a local storage for each type Use
59. itutive Model c Input c eps Current strains at point c theta Trace of strain at point c td Temperature change c d Program material parameters c ud User material parameters c hn nh History terms at point t_n c hi nh History terms at point t_n 1 c nh Number of history terms c ii Current point number c istrt Start state O elastic 1 last solution c isw Solution option from element c Output c sig 6 Stresses at point c dd 6 6 Current material tangent moduli do B use onun qora on GS a e ee implicit none integer nh ii istrt isw i realx8 td 1 8 eps theta d ud hn nh hi nh sig 6 dd 6 6 c Dummy model sig ud 1 eps if isw eq 14 the Set initial values for history parameters None needed for this model Compute tangent and stress else do i 1 6 dd i i ud 1 sig i ud 1 eps i end do endif end Figure 4 4 Sample UMATLn module for small deformation CHAPTER 4 USER FUNCTIONS 25 reference and current coordinates for the constitutive point are passed in common block elcoor which contains the values in realx8 xref XCUr common elcoor xref 3 xcur 3 Internal variable storage and use For constitutive equations with additional internal variables that evolve in time users must define entries for the h1 array The number of entries available in the array for each evaluation i e each quadrature point i
60. l data values to input The value for parameter nt or nn may be zero Thus the use of errck tinput text 0 td nn CHAPTER 2 DATA INPUT AND OUTPUT 7 is equivalent to errck pinput td nn Text variables may be converted to numerical REAL 8 form using the subroutine call call setval text nc td where text is a string with nc characters and td a REAL 8 variable The text string can contain any parameters expressions or numerical constants which evaluate to a single value 2 2 Array Outputs Two subprograms exist to output arrays of integer and real double precision data The routine MPRINT is used to output real data and is accessed by the statement call mprint array nrow ncol ndim label where array is the name of the array to print nrow and ncol are the number of rows and columns to output ndim is the first dimension on the array and label is a character label which is added to the output For example the statements real 8 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 exc
61. lacement in direction zi is an acceleration in direction x v b is a loading in direction x per unit length and des is a virtual strain along the truss axis CHAPTER 5 ADDING ELEMENTS 49 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 K Lais t i 1 where t is time u is the displacement vector with components u l is a unit vector along the axis of the member with direction cosines l defined by On 22 4 Os L E gt xs ra i 1 and i Z are the coordinates of nodes 1 and 2 respectively The displacement components are interpolated on the 2 node truss member as i s 1 un ult in which uj uo 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 g Ou Ess i OS os Using the interpolations for the displacement components yields 1 d Ess T2 Ag 1 1 where Az Z Til 1 and Au Uj2 Ui Thus in matrix form the strain is 1 12 4 Ess a Az Aa d Wig 1 CHAPTER 5 ADDING ELEMENTS 50 Using the above displacement interpolations the variational equation for the truss may be expressed in matrix form as ent un Buso 7 4 ee oss Ads m Al elas Pad eM ES a R
62. m 1 dim 2 dim 3 Description CMASn 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 arrays IE and IEDOF define additional information required to process each element Tables 3 4 3 5 and 3 6 describe the use of individual entries in the arrays IX IE and IEDOF 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 setvar palloc 111 TEMP1 palloc 112 2 2 1 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 CHAPTER 8 ALLOCATING ARRAYS NAME Num dim 1 dim 2 dim 3 Description ANGL 46 nen Angle LD 35 nst Assembly nos P 35 nst Element vector S 36 nst nst Element matrix TL 39 nen Te
63. mmand line input parameter name c uprt Flag Output results if true c Outputs c none Users responsible for outputs to arrays etc c implicit none include umaci h 1 Contains UCT variable character tx 15 logical uprt Set name mesi to user defined if pcomp uct mesi 4 then uct Set user defined command name else User execution function statements follow end if end Figure 4 1 Sample UMESHn module CHAPTER 4 USER FUNCTIONS 17 The ability to get array names as shown in Chapter 3 can be used to develop user routines for input of coordinates element connections etc With this facility it is possible to develop an ability to directly input data prepared by other programs which may be in a format which is not compatible with the requirements of standard FEAP mesh commands 4 1 1 Command line TX data It is possible to include up to 8 data items on the command line for user functions All the data is passed to the UMESHn functions by the character array TX 15 and may be used to control actions in the function If the information is of type character it may be used directly however if it is numeric it must be converted within the UMESHn 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 UME
64. mperature UL 41 ndf nen 6 Solution array XL 44 ndm nen Coordinates Table 3 3 Element Array Names Numbers and Sizes NAME Description IX 1 ei Global node 1 Rita to IX nen e Global node nen IX nent1 e IX nen 2 e IX nen 3 e IX nen 4 e IX nen 5 e IX nen 6 e IX nent7 e H1 history data pointer H2 history data pointer H3 history data pointer Lagrange multiplier tag Lagrange multiplier data pointer Time integrator 0 implicit gt 0 explicit Element type FE lt 0 IGA gt 0 IX nen1 ei IX nen1 1 e IX neni 2 e Element material type number Element region number default 0 Active region gt 0 Inactive region lt 0 Active deactive start Table 3 4 Element connection array IX use for element e 12 CHAPTER 8 ALLOCATING ARRAYS IE nie ma IE nie 1 ma IE nie 2 ma IE nie 3 ma IE nie 4 ma IE nie 5 ma IE nie 6 ma IE nie 7 ma IE nie 8 ma IE nie 9 ma IE nie 10 ma NAME Description IE 1 ma Plot shape dimension 0 1 2 3 O no plot 1 line 2 surface 3 solid IE 2 ma Rigid material number Number history variables element NH1 and NH2 Element material type number ELMTO1 1 etc Element material type identifier 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 fo
65. ng model energy parameter 176 ag Fung model energy parameter 177 az Fung model energy parameter 178 ag Fung model energy parameter 179 ag Fung model energy parameter 180 181 Viscoplastic rate parameters 182 Nodal quadrature parameters 183 Ba Mz Mg mass scaling factor 184 ic Estimate on maximum wave speed 185 Augmentation switch lt on off gt 186 Augmentation explicit indicator 187 Implicit 0 Explicit 1 element integration Table 5 8 Material Parameters CHAPTER 5 ADDING ELEMENTS 43 Parameter Name Description 188 Number stress components in rod elements 189 Nurbs flag 190 192 Nurbs quadrature values direction 193 tmat Thermal material numbers 194 ietype Element type 195 T frac Fraction of work to heat 196 q prop Proportional load factor for pressure loading 197 198 Body patch loading values 199 Axisymmetric 1 d Plane stress in thickness 200 nsiz Size of modulus or compliance array 201 236 Anisotropic Modulus or Compliance array 237 239 Unused 240 0 Element based 1 nodal based formulation 241 Number of active element degrees of freedom 242 250 Unused Table 5 9 Material Parameters and Viltn At 1 y ailtn vai tas with u v and a being the set of displacements velocities and accelerations at node i respectively A Newton method is commonly adopted to solve a non l
66. nt axial force eps I Element axial strain could be placed anywhere after the stress and strain are defined These values would be output by using a solution command sequence such as batch tplot end stress nn i saves force for element nn stress nn 2 saves strain for element nn show writes tplot items to output file 5 5 Projection of element variables to nodes The STREss NODE solution command and the PLOT STRE command require a projection of element variables to nodes A continuous stress field is assumed to obtain the nodal values Accordingly o where o 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 ca nodal values of the projected quantity The projection routine uses a diagonal weight matrix to project the values For simple elements the matrix is computed by a procedure identical to mass lumping For example Man d Na dO Q defines a row sum form of projection 3 Using the above results in the set of equa tions and a least square fit with the finite element values 6 gives the equation set Moa Ca a Naa dt 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 ject
67. of time steps niter Number of iterations current step naugm Number of augments current step titer Total iterations taubm Total augments iaugm Augmenting counter iform Number residuals in line search dm Element proportional load n Current element number ma Current element material set mct Print counter iel User element number nel Number nodes on current element tt Element stress values for TPLOt bpr Principal stretch ctan Element multipliers ut Element user values for TPLOt Table 5 3 FEAP common block definitions 5 2 Non linear Transient Solution Forms Before describing the steps in developing an element we summarize first the basic 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 equations expressed in terms of nodal displacements velocities and accelerations given by u t t and u t respectively We denote the differential equation for node i as the CHAPTER 5 ADDING ELEMENTS 38 Variable Definition nh1 Pointer to t history data nh2 Pointer to t 1 history data nh3 Pointer to element history ior Current input logical
68. outines described below 2 Node Element 3 Node Element Figure 6 1 Line type elements in FEAP library Lagrangian interpolation in one dimensional isoparametric forms may be obtained us ing the call CALL SHP D 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 1 the shape function N In calculations integrals are represented as ANo Ng ds JING NG OX IAC de 6 5 and quadrature may be used for evaluation Calculation of natural coordinate derivatives only may be obtained with the call CHAPTER 6 UTILITY ROUTINES 75 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 NY n NY z N 0i N 0 6 6 is obtained using the call CALL SHP1DH S LEN SHPW SHPT where Parameter Description 5 Natural coordinate LEN Length of the element 2 node SHPW 4 2 Shape func
69. p lct off 3 then solver true Sets flag for FEAP solvers any other statements needed else solver false Sets flag for user solver any other statements needed endif endif end and the module usolve f 87 CHAPTER 7 ADDING A USER SOLVER subroutine usolve flags b 88 00 0 Purpose Solver interface for SuperLU Inputs flags 1 Allocation and or initialization phase flags 2 Perform factorization for direct solutions flags 3 Coefficient array unsymmetric flags 4 Solve equations flags 5 Purge storage of pointers b RHS vector Outputs flags 5 True if error occurs for factor solve only Ee r a implicit none logical flags real 8 b Presolve setups if flags 1 then Solution steps for assembled equations else Factor equations if flags 2 then endif Perform solve if flags 4 then endif Purge storage in factor if flags 5 then endif endif end Bibliography 1 G G Weber A M Lush A Zavaliangos and L Anand An objective time integration procedure for isotropic rate independent and rate dependent elastic plastic constitutive equations International Journal of Plasticity 6 701 744 1990 21 O C Zienkiewicz and R L Taylor The Finite Element Method volume 1 McGraw Hill London 4 edition 1989 3 O C Zienkiewicz and R L Taylor The Finite Element Met
70. placement are given by d unj s t mn als do Donguls t i 1 and the derivative by r gt ng a Collecting terms and combining with previously defined quantities the virtual strain may be written as Oou o 18 song After differentiation of the displacement field the discrete form of the virtual strain is de where given by 1 de uu du i Substituting the above virtual strain expression into the weak form gives the modified residual expression 4 1 Ji L 2 1 l n nal a fa saz ll s 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 g B and O u du 2 n n n n SS If the n are constructed as column vectors then the tensor product becomes a matrix defined as d c T G n n np n n n n With these definitions the tangent matrix for the non linear problem is given as EA di Oss Gij Gij Ki 9 9 s l G Gl 72 CHAPTER 5 ADDING ELEMENTS 63 Notice that for the linear problem Az L Ji thus the only difference betvveen 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 is sometimes called the geometric stiffness part CHAPTER 5 ADDING EL
71. r PROTxx or UROT xx Get tangent from element if 0 if gt 0 numerically differ entiate residual to obtain tangent Equation number for element Lagrange multiplier Partition number for element Lagrange multiplier Global equation number Table 3 5 Element control array IE use for material number ma NAME Description IEDOF 1 i ma IEDOF ndf i ma Degree of freedom 1 for node i of material ma to Degree of freedom ndf for node i of material 13 Table 3 6 Element degree of freedom assignment array IEDOF use for material number ma CHAPTER 8 ALLOCATING ARRAYS 14 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 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 tru
72. r problems in solid mechanics the linear momenta are stored as follows The linear momenta are computed as p f eve Q n l x x pda Q the angular momenta as the kinetic energy K dQ Q CHAPTER 5 ADDING ELEMENTS 61 the stored energy as U W C dQ Q and the work by external loads as V eH X Fea ar r The value of the displacement and velocity at the current time t are passed in ul i j 1 and ul i j 4 respectively Note that this is true no matter which time integration algorithm is specified 5 9 A Non linear Theory for a Truss A simple non linear theory for a two or three dimensional truss which may undergo large displacements for which the strains remain small may be developed by defining the axial strain approximation in each member as E Na 1 r a Ess Os 2 j 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 OUn gt Os GD An algorithm to define the two orthogonal unit vectors which are normal to the member may be constructed by taking V EL where k is a direction for which a minimum value of the direction cosine l exists for a 2 dimensional problem defined in the z z plane v may be taken as e Now vx1 xd Di and no Ixn CHAPTER 5 ADDING ELEMENTS 62 Using these vectors the two normal components of the dis
73. ress strain etc are developed from task 8 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 Unta j 2 Increment ul Un j 3 Increment 4 j 5 Rate all 6 Rate un 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 39 It is not necessary to implement all other tasks in an element however for those tasks that are not implemented it is important that the element routine not perform any calculations Thus if the form of the branch is programmed as an IF THEN ELSE construct as shown in Fig 5 1 then the ELSE should not carry out any operations unless all options for ISW are programmed Similarly if the element is programmed using a SELECT CASE form as isw Task Type Description Access Command 0 O Output label
74. routine as S ctan D 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 CHAPTER 5 ADDING ELEMENTS 47 Parameter Description ctan 1 cj Multiplier of s matrix for ul i j 1 terms e g stiffness matrix multiplier ctan 2 co 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 10 Tangent Parameters 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 isw we consider next the problem of a 2 node linear elastic truss element for small deformation applications The element is described in sufficient generality to permit solution of both two and three dimensional truss problems 5 3 1 Theory for a
75. rray 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 912 Sis 5 S r 922 S23 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 is passed to the subprogram as ndf Thus SS 5G mala in which is an array coefficient for nodal pair i j for the degree of freedom pair a b m 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 Ro R Rs CHAPTER 5 ADDING ELEMENTS 46 where the entries in each submatrix are given as i Ri 2 i 2 R R3 RA The two dimensional form places the entries R as columns Accordingly R R R Rs SC The two forms for defining the residual r are equivalent based on the Fortran ordering of
76. rs 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 locations The element history data associated with tn starts at the memory address of the pointer for NH1 using the double precision dummy array HR in blank common similarly data for t 1 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 Parameters of type NH3 may not be used in material model routines UMATLn 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 CHAPTER 5 ADDING ELEMENTS 56 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 nh1 nh3 6 10 reserves 6 words of NH1 and NH2 data and 10 words of NH3 data for each element with the current material number Care should be taken to minimize the number of history variabl
77. s nh The value for nh is defined by the parameter n1 in module UMATIn see Fig 4 3 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 array are copied to the hn x array each time the command statement TIME is issued in a solution 4 2 3 Auto time step control The solution command AUTO MATErial rvalu 1 rvalu 2 rvalu 3 initiates an attempt to control the solution process by a variable time stepping algo rithm based on a user set value in the material constitution The value to be set is named rmeas which is passed between constitution and solution modules in the labeled common real 8 rmeas rvalu CHAPTER 4 USER FUNCTIONS 26 subro
78. sidual 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 5 4 Additional Options in Elements FEAP permits some additional options to be included within element tasks 5 4 1 Task 1 Options Setting active equations Often it is necessary to use several element types to perform an analysis For example it may be necessary to use both truss and frame bending resistant elements to perform an analysis As developed in Section 5 3 the truss element has one degree of freedom for each spatial dimension whereas the frame element must have additional unknowns to represent the bending behavior For nodes connected only to truss elements it is not necessary to have the additional degrees of freedoms 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 the degrees of freedom with values greater than ndm the spatial dimension of the mesh CHAPTER 5 ADDING ELEMENTS 52 do i ndmti ndf
79. ss of an array in a calling statement to also be 64 bits in length For example use of integer ioff ioff np 111 numnp call submat hr ioff would cause an error since the pointer ioff is only 32 bits in length To avoid this problem it is necessary to either declare ioff to be 64 bits long as integer 8 ioff or use one of the FEAP include files p int h defining the integer type array fp 10 or p point h defining the integer type scalar point A type of correct length is controlled at compile time by which subdirectory is used i e integer4 or integer8 Using this scheme permits direct reference to either real 8 or integer arrays in pro gram modules without need to pass arrays through arguments of subprograms A subprogram PALLOC controls the allocation of all standard arrays in FEAP defined by the np pointers and a subprogram UALLOC permits users to add allocation for their own arrays defined by the pointers up The basic use of the routines is provided by an instruction setvar palloc number NAME length precision or setvar ualloc number NAME length precision where setvar palloc and ualloc are logical types number is an integer number of the array NAME is a 5 character name of the array length is the number of words of storage needed for the array and precision is the type of array to allocate 1 for integer and 2 for real 8 types Upon initial assignment of any array its values are set
80. ta is normally provided as numerical data without use of any parameters or expressions CHAPTER 1 INTRODUCTION 2 the input time may be reduced by setting the value of the logical variable COFLG in feap84 f to false FEAP will automatically switch to parsing mode if any record contains non numerical data item It is also possible to use the PARSe and NOPArse commands to set the appropriate mode of data input In Windows versions it is sometimes desirable to obtain the input file name from a pop up menu This is accomplished by setting the parameter CIFLG to true 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 feap84 f 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 The last parameter which may be set in the feap84 f module is the level for displaying
81. ther than giving equivalent common statements directly If later releases of the FEAP program revise contents in a common block it will only be necessary to recompile the user routine rather than change all the common statement definitions Also by defining the correct path in the makefile in on your compiler it is not necessary to modify any routine when switching from 32 bit machines to 64 bit machines Chapter 2 DATA INPUT AND OUTPUT FEAP includes utilities to perform input and to output small arrays of data Users are strongly encouraged to use the input utilities but often may wish to use their own utilities to output data 2 1 Parameters and Expressions The subroutines PINPUT and TINPUT are input subprograms used by FEAP to input each data record They permit the data to be in a free form format with up to 255 characters on each record as well as to employ expressions parameters and numerical representations for each data item These routines also should be used to input data in any new program module developed The PINPUT routine returns data to the calling subprogram in a double precision array The following statements may be included as part of the routine performing the input subroutine xxx logical errck pinput include iofile h ior iov ilg unit numbers real 8 td 5 1 if ior lt 0 write 3000 errck pinput td 5 if errck go to 1 The parameters defined in the include file common block are C
82. tions 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 xs 3 SHPW 3 i SHPT 3 i are third derivatives e g Nj srx and 4 SHPW 4 i SHPT 4 i are shape functions e g N 6 2 2 Shape functions in two dimensions The shape funtions for two dimensional triangular elements as shown in Fig 6 2 and quadrilateral elements as shown in Fig 6 3 may be computed using the shape function routines described below CHAPTER 6 UTILITY ROUTINES 76 Two dimensional Cy isoparametric interpolation on quadrilaterals of linear quadratic and cubic order may be obtained using the subprogram call CALL SHP2D SS XL SHP XJAC NDM NEL IX FLG where Parameter Description 56 2 Natural coordinates 7 XL NDM NEL Element coordinates in local order NDM Spatial dimension mesh 2 or 3 NEL Number nodes on element 4 9 12 16 IX NEL Element global node numbers FLG Return 7 derivatives if true or x y derivatives if false SHP 3 NEL Shape functions and derivatives XJAC Jacobian transformation from z y to n The array SHP stores the values in the order SHP 1 i derivative with respect to or x SHP 2 i derivative with respect to 7 or y SHP 3 i shape function Two dimensional Co isoparametric interpolation on triangles of linear
83. tory Variables 5 6 1 Assigning amount of storage for each element 5 6 2 Accessing history data for each element 15 15 17 18 20 20 25 27 29 29 CONTENTS 5 7 Elements with Finite Rotation Parameters 02 5 7 1 Nodal rotation treatment UROTmm subprogram 5 7 2 Local nodal rotation treatment 200 5 8 Energy Computation 5 9 A Non linear Theory for a Truss ss cm o 2 x hoe ee sil eS 6 Utility routines 6 1 Numerical quadrature routines se ps pr ep Epa 611 One dimensional quadrature pus oz WA 6 1 2 Two dimensional quadrature bk eee bat Ds 6 1 3 Three dimensional quadrature 6 2 Shape function subprograms cu ane dm a a a eS B 6 2 1 Shape functions in one dimension 6 2 2 Shape functions in two dimensions 6 2 3 Shape functions in three dimensions 6 3 Eigenvalues for 3 x 3 matrix ha Dada m x oa AS d R 6 4 Plot routines 6 4 1 Mesh plots 6 4 2 Element data plots iis 4 gs sala 4 gs kee arad a oer A ASG 6 4 3 Other user plots 7 Adding a user solver 57 58 59 60 61 70 70 70 71 72 73 74 75 78 80 81 81 83 89 87 List of Figures 4 1 4 2 4 3 4 4 4 5 4 6 4 7 4 8 5 1 5 1 5 2 5 3 5 4 5 5 5 6 5 7 5 8 6 1 6 2 6 3 6 4 6 5 Sample MESHTm dl le k
84. unit iow Current output logical unit nph Pointer to global projection arrays ner Pointer to global error indicator erav Element error value j int J integral values ndf Maximum dof node ndm Mesh space dimension nent Dimension 1 on IX array nst Size of element matrix nneq Total dof in problem ttim Current time dt Current time increment ci Integration parameters hr Real array data mr Integer array data Table 5 4 FEAP common block definitions residual equation To solve for the nodal displacements velocities and accelerations it is necessary to introduce an algorithm to integrate the nodal quantities in time specify a constitutive relation and develop an algorithm to solve a possibly non linear problem In FEAP the integration method for nodal quantities is taken as a one step algorithm with each quantity defined only at discrete times t Accordingly we have displace ments u tn with velocities and accelerations denoted as tn So Vi tn and tn ai tn A typical example for an integration algorithm for these discrete quantities is New mark s method where 1 uti Uj tny Atvi try AP IS B a try Bai tr4iy CHAPTER 5 ADDING ELEMENTS Parameter Name Description 11 Young s modulus A IR Poisson ratio ala Thermal expansion coefficient 419 Mass density 51 Quadrature order for arrays 6 Quadrature order for outputs Tia Mass interpolat
85. utine umatl1 f detf td d ud hn h1 nh ii istrt sig dd isu AA NAA ne AA Sa eee Purpose User Constitutive Model Input 3 3 Deformation gradient finite deformation detf Determinant of deformation gradient td Temperature change d Program material parameters ud User material parameters hn nh History terms at point t_n hi nh History terms at point t n 1 nh Number of history terms ii Current point number istrt Start state O elastic 1 last solution isw Solution option from element Output sig 6 Stresses at point dd 6 6 Current material tangent moduli eee q er ic m nm 1 implicit none integer nh istrt ii isu realx8 td realx8 3 3 detf x d k ud hn nh h1 nh 1 8 sig 6 dd 6 6 Model if isw eq 14 then 1 Set any initial values for history else Compute model tangent and stress endif end Figure 4 5 Sample UMATLn module for finite deformation CHAPTER 4 USER FUNCTIONS 27 logical aratfl common elauto rmeas rvalu 3 aratfl The three parameters may be used in defining an acceptable value for rmeas The algorithm coded monitors the solution during a standard iteration process set by for example LOOP n TANG 1 NEXT If during any iteration up to n the value of rmeas exceeds a value of 2 rmeas 0 at the start of the loop a new value of At is immediately set to Atnew 0 85 At rmeds and the iteration process is
86. ws 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 CHAPTER 4 USER FUNCTIONS 31 Number Color Black White Red Green Blue Yellow Cyan Magenta Orange Coral 10 Green Yellow 11 Wheat 12 Royal Blue 13 Purple 14 Aquamarine 15 Violet Red 16 Dark Slate Blue 17 Gray 18 Light Gray QO ND ot Ah ON HO Table 4 1 Color Table for Plots 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 5 ADDING ELEMENTS FEAP permits users to add their own element modules to the program by writing a single subprogram called subroutine elmtnn d ul xl ix tl s r ndf ndm nst isw where nn may have values between 01 and 50 Each element subprogram must be added before loading the FEAP library since dummy subprograms are included in the library to avoid unsatisfied externals The basic structure for an element routine is shown in Figures 5 1 and 5 1 Information is provided to the element subprogram through data passed as arguments and data passed in common blocks The data passed as arguments consists of eleven subroutine elmtnn d ul xl ix tl s r ndf ndm nst isw
87. ya ao m ak Bl Q l g d k l 16 Sample UMESH module 4953 45 ed ma ii 19 Sample UMATI1 module ms Kna Sele e xd S y E SS 21 Sample UMATLn module for small deformation 24 Sample UMATLn module for finite deformation 26 Sample UMANLn module 28 Sample UMACRn module et ae ed eee ee Gree KAT KAG 29 Sample UPLUTn module e e Bla AA tele e i 30 FEAP Element Subprogram Part 1 32 FEAP Element Subprogram Part 2 ga a aa Bea ng 33 FEAP Element Common Blocks sales eke E AA 64 FEAP Element Common Blocks using Includes 65 2 Node Truss Element s a d RU L kee EE 65 Element residual for two node truss 66 Truss Tangent Matrix Part sl x gata ack de Pooh ea Soa eg 67 Truss Tangent Matrix Part 2 E A eo eR GA 68 Element variable projection routine 69 Line type elements in FEAP library 74 Triangular surface type elements in FEAP library 76 Quadrilateral surface type elements in FEAP library 77 Tetrahedron solid type elements in FEAP library 78 Brick solid type elements in FEAP library 79 iii List of Tables 3 1 3 2 3 3 3 4 3 5 3 6 4 1 5 1 5 2 5 3 5 4 5 9 5 6 5 7 5 8 5 9 5 10 5 11 5 12 6 1 6 2 6 3 6 4 Mesh Array Names Numbers and Sizes
88. yield stress Mises 43 8 Exponential hardening rate 44 Hiso Isotropic hardening modulus linear 45 Hkin Kinematic hardening modulus linear 46 Yield flag AT Bi Orthotropic thermal stress 48 z Orthotropic thermal stress 49 B3 Orthotropic thermal stress 50 Error estimator parameter 51 yy Viscoelastic shear parameter 52 Ti Viscoelastic relaxation time 53 Vo Viscoelastic shear parameter 54 T Viscoelastic relaxation time 55 13 Viscoelastic shear parameter 56 73 Viscoelastic relaxation time 57 nvis Number of viscoelastic terms 1 3 58 Damage limit 59 Damage rate 60 k Penalty parameter Table 5 6 Material Parameters 40 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 Tl gi Ground acceleration factor 72 92 Ground acceleration factor 73 93 Ground acceleration factor 74 py 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 TS ay Rayleigh damping stiffness ratio 79 Plate Shell Rod shear activation flag 80 Method

FEAP Programmer Manual - Department of : Civil and

Contents

Download Pdf Manuals

Related Search

Related Contents