Home

Calfem - Chair of Applied Mathematics / Numerical Analysis

1. Table 1 Analogous quantities 5 2 1 ELEMENT Interpretations of the spring element Problem type Quantities Designations gk k spring stiffness spring 5 yyy Saplacemens Pring P element force AMi SPs L length u P u P E modulus of elasticity B A area of cross section ar E A L f u displacement N N P element force 4 N normal force L length Qo WT Q TQ thermal conductivity Thermal e r conduction r eee eae Q element heat flow L f Q internal heat flow R R resistance Electrical U a U U potential circuit T b I element current T I internal current L length Ground QO Q gt QO k permeability o e water 7 piezometric head flow xe Q element water flow L Q internal water flow L length Pipe D j i p f Di i pipe diameter network u viscosity uwx _Q laminar O Q p pressure flow 7 l Q element fluid flow Q internal fluid flow Table 2 Quantities used in different types of problems ELEMENT 5 2 2 The following functions are available for the spring element Spring functions springle spring1s Compute element matrix Compute spring force 52 3 ELEMENT springle Spring element Purpose Compute element stiffness matrix for a spring ele
2. setts bday 7 o4 A 4 24 24 4 24 24 4 A4 4 where 3 2 and Ly y3 yt ELEMENT 57 4 Plate element platrs Purpose Compute section forces in a rectangular plate element Syntax es et platrs ex ey ep D ed Description platrs computes the section forces es and the curvature matrix et in a rectangular plate element The section forces and the curvatures are computed at the center of the element The input variables ex ey ep and D are defined in platre The vector ed contains the nodal displacements a of the element and is obtained by the function extract as ed af u uz w2 The output variables es MT v Msa Myy Mery Vez Vue et k Kaz Kyy Kay contain the section forces and curvatures in global directions 5 7 5 ELEMENT platrs Plate element Theory The curvatures and the section forces are computed according to k ky BG a Kay Mis M My De Moy Voz SE vel o where the matrices D B C and af are described in platre and where o o ony E cH Ox Oy a Oy Ox ELEMENT 5 7 6 6 System functions 6 1 Introduction The group of system functions comprises functions for the setting up solving and elimi nation of systems of equations The functions are separated in two groups Static system functions Dynamic system functions Static system functions concern the li
3. B DBdAC A C717 f q dA A where the constitutive matrix 3 D D 12 and where D is defined by D The evaluation of the integrals for the rectangular plate element is based on the displacement approximation w x y and is expressed in terms of the nodal variables Uj U2 U12 AS w x y N a N C a where Nafi Y x xy y x xy 1 a BR a ab bB a a b ab 0 0 1 O a 2b 0 a 2ab 0 1 0 2a b 0 3a 2ab 1 a Sb a Sab b a a b ab 0 0 1 0 a 2b 0 a 2ab cC 0 1 0 2a b 0 3a 2ab b A a b ab DP a ab ab 0 0 1 0 a 2b 0 a 2ab 0 1 0 2a b 0 8a 2ab b 1 a b a ab b a a a b ab 0 0 1 0 a 2b 0 a 2ab 0 1 0 2a b 0 3a 2ab b e E af u U a u2 and where 1 1 a 5 ts 1 and b 5 Ys ya The matrix B is obtained as 000200 6 4 0 0 bry B VvVN 000002 0 0 2z y 0 000020 0 4r 4y 0 6x 5 7 3 ry y xy ry b a b ab 3b a 3ab 0 3a b 8 b a b ab 3b a Sab 0 3a2b b b ab ab 3b a Bab 0 3a b b b a b ab 3b lt a 3ab 0 3a b b 0 6xy 6y ELEMENT platre Plate element where 32 Ox 32 dy a OxOy lt 2 Evaluation of the integrals for the rectangular plate element is done analytically Computation of the integrals for the element load vector ff yields At ig ge Ws lige S a SL 3 S a fe 9 1
4. Ey 2 ep1 N 2 Ke3 beam2 Ex 3 Ey 3 ep3 N 3 EXAMPLES 9 5 6 Nonlinear analysis exn2 K assem Edof 1 K Ke1 K assem Edof 2 K Ke2 K assem Edof 3 K Ke3 bc 1 0 2 0 3 0 10 0 11 0 12 O a Q solveq K f bc Ed extract Edof a esi beam2gs Ex 1 Ey 1 ep1 Ed 1 N 1 es2 beam2gs Ex 2 Ey 2 ep1 Ed 2 N 2 es3 beam2gs Ex 3 Ey 3 ep3 Ed 3 N 3 NO N N esi 1 1 es2 1 1 es3 1 1 if n gt 20 disp The solution doesn t converge return end end Kred red K bc 1 if det Kred lt 0 disp Determinant lt 0 buckling load passed break end disp Alpha num2str alpha is OK int2str n iterations are performed disp deform j a 4 M j Q 3 loadfact j alpha bmode a end 9 5 7 EXAMPLES exn2 Nonlinear analysis The following text strings are produced by the m file Alpha 1 is OK 3 iterations are performed Alpha 1 1 is OK 3 iterations are performed Alpha 1 2 is OK 3 iterations are performed Alpha 1 3 is OK 3 iterations are performed Alpha 2 4 is OK 3 iterations are performed Alpha 2 5 is OK 3 iterations are performed Alpha 2 6 is OK 3 iterations are performed Alpha 2 7 is OK 3 iterations are performed Alpha 2 8 is OK 3 iterations are performed Alpha 2 9 is OK 3 iterations are performed Alpha 3 is OK 3 iterations are performed Alpha 3 1 is OK
5. 13 ELEMENT beam2ws Two dimensional beam element Theory The section forces at the ends of the beam are obtained from the element force vector P N V M N V M computed according to P K K Ga The matrices K and G are described in beam2e and the matrix K is described in beam2w The nodal displacements T a u U2 UZ U4 U5 ug are shown in beam2e Note that the transpose of a is stored in ed ELEMENT 5 6 14 Two dimensional beam element beam2g Purpose Compute element stiffness matrix for a two dimensional nonlinear beam element X29 Xpy Syntax Ke beam2g ex ey ep N Ke fe beam2g ex ey ep N eq Description beam2g provides the global element stiffness matrix Ke for a two dimensional beam element with respect to geometrical nonlinearity The input variables ex ey and ep are described in beam2e and N N contains the value of the predefined normal force N which is positive in tension The element load vector fe can also be computed if a uniformly distributed transverse load is applied to the element The optional input variable eq qz then contains the distributed transverse load per unit length qx Note that eq is a scalar and not a vector as in beam2e 5 6 15 ELEMENT beam2g Two dimensional beam element Theory The element stiffness matrix K stored in the variable Ke is computed according to K G K G where G is
6. 3 4 ep ptype t n ey y y2 ys ya p typ Dy Di2 Dig Diu Dis Di6 D Dog Do3 D D D D D D 7 11 12 13 B D3 D32 D33 D34 D35 D36 D Da Dog D 3 or D Day Dap Dag D44 Das D46 D3 D32 D33 Ds1 Ds2 Ds3 Dsa Dss Dose Dei Dez Des Dea Des Deo If different D matrices are used in the Gauss points these D matrices are stored in a global vector D For numbering of the Gauss points see eci in plani4s D pel Da ELEMENT 5 5 20 Two dimensional solid elements plani4e If uniformly distributed loads are applied to the element the element load vector fe is computed The input variable by containing loads per unit volume by and b is then given Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT DB tdA A f NTbtdA A with the constitutive matrix D defined by D and the body force vector b defined by eq The evaluation of the integrals for the isoparametric 4 node element is based on a displacement approximation u 7 expressed in a local coordinates system in terms of the nodal variables u1 uo ug as ulg n N af where Ui E e NE 0 NE 0 NE 0 NE 0 o ge vel A 0 A ON aN TS Us The element shape functions are given by e 1 e 1 Ny 70 d m Ny z l 1 7 1 1 Ny 70 8 G 9 Ne AE The matrix B is obtained as o Da 0 2 S o
7. dt 0 002 T 1 boundary condition initial condition G 0 0 0 1 0 02 0 2 0 01 0 3 0 0 T 0 t g gfunc G dt bc zeros 4 1 length g bc 1 1 g bc 2 1 2 bc 3 1 3 bc 4 1 14 d0 zeros 15 1 v0O zeros 15 1 output parameters ntimes 0 1 0 1 1 nhist 1 4 11 time integration parameters ip dt T 0 25 0 5 10 3 ntimes nhist time integration k sparse K m sparse M Dsnap D V A step2 k m d0 v0 ip bc 4 plot time history for two DOF s figure 1 plot t D 1 t D 2 t D 3 grid xlabel time sec ylabel displacement m title Displacement time at the ist 4th and 11th gt degree of freedom text 0 2 0 022 solid line bottom vertical beam 9 4 11 EXAMPLES Dynamic analysis gt x direction text 0 2 0 017 dashed line center vertical beam gt x direction text 0 2 0 012 dashed dotted line center gt horizontal beam y direction plot displacement for some time increments figure 2 clf axis equal hold on axis off magnfac 20 title Snapshots sec magnification 20 for i 1 5 Ext Ex i 1 3 eldraw2 Ext Ey 2 3 0 Edb extract Edof Dsnap i eldisp2 Ext Ey Edb 1 2
8. gt gt a ans 5 12 3 Addition and subtraction of matrices gt gt A B ans 3 11 T 9 10 4 gt gt A B ans 1 1 1 5 6 4 Note that if the result of an operation is not assigned to a specific variable the answer is temporarily stored in the variable ans 9 2 5 EXAMPLES exi2 MATLAB introduction Multiplication of matrices gt gt a b ans 17 gt gt a b ans 5 0 20 12 0 48 3 0 12 gt gt A B ans 16 9 4 68 46 24 34 23 12 gt gt C B A C 44 60 19 30 To perform arithmetic operations matrix dimensions must agree gt gt D A B Error using gt Inner matrix dimensions must agree The inverse of a square matrix gt gt inv C ans 0 1667 0 3333 0 1056 0 2444 EXAMPLES 9 2 6 MATLAB introduction exi2 The determinant of a square matrix gt gt det C ans 180 An array or element by element arithmetic operation is denoted by a period pre ceding an operator Examples are element by element multiplication division and powers gt gt a b Matrices in element by element operations must obviously have the same dimensions Mathematical functions applied to arrays are also evaluated element by element gt gt sin 0 1 2 3 4 pi 2 ans 92T EXAMPLES exi3 MATLAB introduction Purpose Show how to handle script files and function files Description When starting a MATLAB session the
9. is an array with the same number of rows as the length of x and y text marks each point with the corresponding row of the string array Note This is a MATLAB built in function For more information about the text function type help text 8 17 GRAPHICS title Purpose Titles for 2D and 3D plots Syntax title text Description title adds the text string text at the top of the current plot Note This is a MATLAB built in function For more information about the title function type help title GRAPHICS 8 18 xlabel ylabel zlabel Purpose x y and z axis labels for 2D and 3D plots Syntax xlabel text ylabel text zlabel text Description xlabel adds text beside the x axis on the current plot ylabel adds text beside the y axis on the current plot zlabel adds text beside the z axis on the current plot Note This is a MATLAB built in function For more information about the functions type help xlabel help ylabel or help zlabel 8 19 GRAPHICS 9 User s Manual examples 9 1 Introduction This set of examples is defined with the ambition to serve as a User s Manual The examples except the introductory ones are written as m files script files and supplied together with the CALFEM functions The User s Manual examples are separated into four groups MATLAB introduction Static analysis Dynamic analysis
10. mie UNIVERSITY CALFEM A finite element toolbox to MATLAB Version 3 3 Department of Mechanics and Materials Structural Mechanics amp Solid Mechanics Department of Mechanics and Materials Structural Mechanics ISRN LUTVDG TVSM 99 9001 SE 1 265 ISSN 0281 6679 CALFEM A finite element toolbox to MATLAB Version 3 3 Copyright 1999 by Structural Mechanics LTH Sweden Printed by JABE Offset Lund Sweden The software described in this document is furnished under a license agreement The software may be used or copied only under terms in the license agreement No part of this manual may be photocopied or reproduced in any form without the prior written consent by the Division of Structural Mechanics CALFEM February 1999 Copyright 1992 99 by the Division of Structural Mechanics and the Department of Solid Mechanics at Lund University All rights reserved CALFEM is the trademark of the Division of Structural Mechanics Lund University MATLAB is the trademark of The MathWorks Inc The Division of Structural Mechanics The Department of Solid Mechanics Lund University Lund University PO Box 118 PO Box 118 S 221 00 Lund S 221 00 Lund SWEDEN SWEDEN Phone 46 46 222 0000 Phone 46 46 222 0000 Fax 46 46 222 4420 Fax 46 46 222 4620 E mail addresses calfem byggmek 1th se CALFEM questions and comments strucmech byggmek 1th se general questions to the departments Homepage http
11. nxn integration points n 1 2 3 are supplied to the function by ep and the thermal conductivities or corresponding quantities kar ky etc are supplied by D ex t1 Z2 3 T4 kez key ey y Yo Y3 ya eitn p f t yx yy If the scalar variable eq is given in the function the element load vector fe is com puted using eq Q where Q is the heat supply per unit volume ELEMENT 5 4 8 Two dimensional heat flow elements flw2i4e Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT DB dA f NT Quad with the constitutive matrix D defined by D The evaluation of the integrals for the isoparametric 4 node element is based on a temperature approximation T 7 expressed in a local coordinates system in terms of the nodal variables T1 75 T3 and T as T 7 Nfae where Ne Nf NE NS N ak T To Ts T The element shape functions are given by Np 71 O n NP FF O 1 4 1 1 Nz 30 1 n N ne 1 n The B matrix is given by o Ch Ein Ox e__ T 1 o e a 2 Oy On where J is the Jacobian matrix ax ar ya 8 On Oy Oy Of On Evaluation of the integrals is done by Gauss integration 5 4 9 ELEMENT flw2i4s Two dimensional heat flow elements Purpose Compute heat flux and temperature gradients in a 4 node isoparametric heat flow element Syntax es et eci flw2i4
12. u uz ug The output variables es o7 Osz Oyy oz Ory ox oy et e See Eyy Ezz Yey Mezl yzl contain the stress and strain components The size of es and et follows the size of D Note that for plane stress 0 and for plane strain o 0 Theory The strains and stresses are computed according to e Bea o De where the matrices D B and a are described in plantce and where the evaluation point x y is chosen to be at the center of the element 5 9 19 ELEMENT plani4e Two dimensional solid elements Purpose Compute element matrices for a 4 node isoparametric element in plane strain or plane stress ls Syntax Ke plani4e ex ey ep D Ke fe plani4e ex ey ep D eq Description plani4e provides an element stiffness matrix Ke and an element load vector fe for a 4 node isoparametric element in plane strain or plane stress The element nodal coordinates 71 y1 2 etc are supplied to the function by ex and ey The type of analysis ptype the element thickness t and the number of Gauss points n are supplied by ep ptype 1 plane stress n x n integration points ptype 2 plane strain n 1 2 3 The material properties are supplied by the constitutive matrix D Any arbitrary D matrix with dimensions from 3 x 3 to 6 x 6 maybe given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 ex 1 we
13. 1 1 1 1 1 1 Orr Oy zz zy Lz yz 2 2 2 2 2 2 es ot Ore Oyy Ozz Ory Otz Oyz n3 n3 n3 n3 n3 n3 LX Tuy Ozz Ory Ozz Oyz 1 1 1 1 1 1 Ers Eyy Ezz Yay Vaz Yyz T Y Z1 2 2 2 2 2 2 E E E T2 Y2 272 et eT rx yy zz Vry Vrz Yyz eci f 3 f 3 i 3 i 3 l Eii Eyy Ezz Yy Yoz Yy Tn Yn Zn contain the stress and strain components and the coordinates of the integration points The index n denotes the number of integration points used within the ele ment cf soli8e 5 9 35 ELEMENT soli8s Three dimensional solid elements Theory The strains and stresses are computed according to e Ba o De where the matrices D B and af are described in soli8e and where the integration points are chosen as evaluation points ELEMENT 5 5 36 Three dimensional solid elements soli8f Purpose Compute internal element force vector in an 8 node isoparametric solid element Syntax ef soli8f ex ey ez ep es Description soli8f computes the internal element forces ef in an 8 node isoparametric solid element The input variables ex ey ez and ep are defined in soli8e and the input variable es is defined in soli8s The output variable f fa fire fiza contains the components of the internal force vector Theory The internal force vector is computed according to f B o av V where the matrices B and o are defined in soli8e and soli8s respectively Evaluation of the integral is done by Gauss integrati
14. 1500 1500 gt gt Ke2 spring1e ep2 Ke2 3000 3000 3000 3000 The element stiffness matrices are assembled into the global stiffness matrix K ac cording to the topology gt gt K assem Edof 1 K Ke2 K 3000 3000 0 3000 3000 0 0 0 0 gt gt K assem Edof 2 K Ke1 K 3000 3000 0 3000 4500 1500 0 1500 1500 9 3 3 EXAMPLES exsl Static analysis gt gt K assem Edof 3 K Ke2 K 3000 3000 0 3000 7500 4500 0 4500 4500 The global system of equations is solved considering the boundary conditions given in bc gt gt bc 1 0 3 0 gt gt a solveq K f bc Element forces are evaluated from the element displacements These are obtained from the global displacements a using the function extract gt gt ed1 extract Edof 1 a edi 0 0 0133 gt gt ed2 extract Edof 2 a ed2 0 0133 0 gt gt ed3 extract Edof 3 a ed3 0 0133 0 EXAMPLES 9 3 4 Static analysis exsl The spring forces are evaluated using the function spring1s gt gt esl spring1s ep2 ed1 esi 40 gt gt es2 spring1s ep1 ed2 es2 20 gt gt es3 spring1s ep2 ed3 es3 40 9 3 5 EXAMPLES exs2 Static analysis Purpose Analysis of one dimensional heat flow Description Consider a wall built up of concrete and thermal insulation The outdoor temperature is 17 C and the temperature inside is 20 C T 20 C surface thermal r
15. 3 9 5 7 C 20 28 46 74 The second type of m files is the function file The first line of these files contains the word function The example below is a function that computes the second and third power of a scalar function b c funcl a 4 function file funci m b axa c axa a 4 end The semi colon prohibits the echo display of the variables on the screen The file can be created using an ordinary editor and it must be saved as funcl m i e the file name without extension must be the same as the function name The second and third power of 2 are calculated by typing gt gt b c funcl 2 producing b 4 Cc 8 See also function and script in Section 7 9 2 9 EXAMPLES exi4 MATLAB introduction Purpose Show different display formats Description Consider the following matrix operation gt gt A 0 4 2 8 gt gt B 3 9 5 7 gt gt C A B 2537 C 0 0079 0 0110 0 0181 0 0292 The result from the computation of C above is shown in the default format displaying four significant decimal digits Other formats can be defined by the command format gt gt format long gt gt C C 0 00788332676389 0 01103665746945 0 01813165155696 0 02916830902641 gt gt format short e gt gt C 7 8833e 03 1 1037e 02 1 8132e 02 2 9168e 02 gt gt format long e gt gt C S 7 883326763894364e 03 1 103665746945211e 02 1
16. Saye aisnast 8 2 bar2e 5 3 2 DAP 0 ceee 5 3 3 bar2s 5 3 5 barge 5 3 6 bar3s 5 3 7 beam2d 5 6 20 beam2ds 5 6 22 beam2e 5 6 2 beam2g 5 6 15 beam2gs 5 6 18 beam2s 5 6 5 beam2t 5 6 7 beam2ts 5 6 9 beam2w 5 6 11 beam2ws 5 6 13 beambe 5 6 24 beam3s 5 6 27 cleat ice scene 2 2 CIP oN omcnpnctneheds 8 3 coordxtr 6 2 3 det cucce 3 7 diat arlenna 3 8 diary eesis 2 3 disp Se Satta See 2 2 4 dmises 4 4 dyna a a 6 3 2 dyna2f wes tae Gund 6 3 3 echo enaa 2 5 CIBC 22 sf act 6 2 5 eldisp2 8 5 eldraw2 8 4 eldia2 8 6 elflux2 8 8 eliso2 8 9 elprinc2 8 10 extract 6 2 6 Mb eeii ari 6 3 4 figure 8 11 TN E 8 12 flw2i4e 5 4 8 flw2i4s 5 4 10 AW2I8E 5 4 11 flw2i8s x2 a2 Srey 5 4 13 flw2qe 162 ces 5 4 6 flw2qs a n 5 4 7 flw2te a n 5 4 3 Wats iaeiae 5 4 5 flw3i8e 5 4 14 flw3i8s 5 4 16 iC eee ee ee eee 7 3 format 2 6 freqresp 6 3 5 Tl nnns as 3 9 gfunc anaana 6 3 6 PUIG inna 8 13 Bel prorcctuh wokus D hold nannan 8 14 hooke cue dante 4 2 Ti cast aid Mosd orinn 7 2 Ty Seenen 6 3 7 insert 6 2 8 y aeii 3 10 length o n 3 11 load oonan 2 8 M
17. a gt gt ed7 extract Edof 7 a gt gt ed9 extract Edof 9 a gt gt Ni bar2s exl eyl ep ed1 Ni 6 2594e 05 gt gt N2 bar2s ex2 ey2 ep ed2 N2 4 2310e 05 gt gt N3 bar2s ex3 ey3 ep ed3 N3 1 7064e 05 gt gt N4 bar2s ex4 ey4 ep ed4 N4 1 2373e 04 gt gt N5 bar2s ex5 ey5 ep ed5 N5 6 9447e 04 gt gt N6 bar2s ex6 ey6 ep ed6 N6 1 7064e 05 gt gt N7 bar2s ex7 ey7 ep ed7 N7 2 7284e 05 EXAMPLES ed2 extract Edof 2 a ed4 extract Fdof 4 a ed6 extract Edof 6 a ed8 extract Edof 8 a edi0 extract Edof 10 a 9 3 14 Static analysis exs4 gt gt N8 bar2s ex8 ey8 ep ed8 N8 2 4132e 05 gt gt N9 bar2s ex9 ey9 ep ed9 N9 3 3953e 05 gt gt N10 bar2s ex10 ey10 ep ed10 N10 3 7105e 05 The largest normal force N 0 62 MN is obtained in element 1 and is equivalent to a normal stress o 250 MPa 9 3 15 EXAMPLES exs5 Static analysis Purpose Analysis of a plane frame Description A frame consists of one horizontal and two vertical beams according to the figure below ee Le ae Ay 45 3e74 m e ane om h 2510e 8 m4 Ao 142 8e 4 m I 33090e78 m E 2108 MPa E P 0 001 MN i go 0 075 MN m The corresponding finite element model consists of three beam elements and twelve degrees of freedom A topology matri
18. bar2e Two dimensional bar element Purpose Compute element stiffness matrix for a two dimensional bar element u4 U3 X22 Syntax Ke bar2e ex ey ep Description bar2e provides the global element stiffness matrix Ke for a two dimensional bar ele ment The input variables ex z 2 ep EA evy v pet supply the element nodal coordinates 71 y1 2 and y2 the modulus of elasticity EF and the cross section area A Theory The element stiffness matrix K stored in Ke is computed according to K G K G where EA 1 1 Nez Tagg 0 0 wl J c 0 O nes ny The transformation matrix G contains the direction cosines T2 T Y2 Y Nez i Nya ja where the length L V 22 1 yo m1 ELEMENT 5 3 2 Two dimensional bar element bar2g Purpose Compute element stiffness matrix for a two dimensional geometric nonlinear bar Ug U3 X Y2 Xpy Syntax Ke bar2g ex ey ep N Description bar2g provides the element stiffness matrix Ke for a two dimensional geometric non linear bar element The input variables ex ey and ep are described in bar2e The input variable N N contains the value of the normal force which is positive in tension Theory The global element stiffness matrix K1 stored in Ke is computed according to K G K G where i lt i 0 0 0 0 _ EA o 0 0 o Nlo 1 0 Se et w ON ele e SO 0 0 0 0 0 1 0 1 Nee
19. n an2 Tn Yn Err Eyy Ezz Yay Vrz Vyz contain the stress and strain components and the coordinates of the integration points The index n denotes the number of integration points used within the ele ment cf plani8e The number of columns in es and et follows the size of D Note that for plane stress 0 and for plane strain o 0 5 9 29 ELEMENT plani8s Two dimensional solid elements Theory The strains and stresses are computed according to e Ba o De where the matrices D B and af are described in plani8e and where the integration points are chosen as evaluation points ELEMENT 5 5 30 Two dimensional solid elements plani8f Purpose Compute internal element force vector in an 8 node isoparametric element in plane strain or plane stress Syntax ef plani8f ex ey ep es Description plani8f computes the internal element forces ef in an 8 node isoparametric element in plane strain or plane stress The input variables ex ey and ep are defined in plani8e and the input variable es is defined in plani8s The output variable ef ff fa fiz fas contains the components of the internal force vector Theory The internal force vector is computed according to f Blo tdA A where the matrices B and o are defined in plani8e and plani8s respectively Evaluation of the integral is done by Gauss integration 5 9 31 ELEMENT soli8e Three dimensional solid element
20. number of timesteps 1 The matrix pbound is organized in the following manner dof time history of the field variable dof time history of the field variable pbound dofm time history of the field variable The dimension of pbound is number of dofs with prescribed field values x number of timesteps 2 The time history functions can be generated using the command gfunc If all the values of the time histories of f or pbound are kept constant these values need to be stated only once In this case the number of columns in f is one and in pbound two It is highly recommended to define f as sparse a MATLAB built in function In most cases only a few degrees of freedom are affected by the exterior load and hence the matrix contains only few non zero entries The computed snapshots are stored in Tsnap one column for each requested snapshot according to ip i e the dimension of Tsnap is number of dofs x nsnap The computed time histories of d and d are stored in D and V respectively one line for each requested degree of freedom according to ip The dimension of D is nhist x number of timesteps 1 6 3 11 SYSTEM step2 Dynamic system functions Purpose Compute the dynamic solution to a set of second order differential equations Syntax Dsnap step2 K C M d0 v0 ip f pdisp Dsnap D V A step2 K C M d0 v0 ip f pdisp Description step2 computes at equal time steps the solution to a second
21. plane stress Syntax es et planrs ex ey ep D ed Description planrs computes the stresses es and the strains et in a rectangular Melosh element in plane strain or plane stress The stress and strain components are computed at the center of the element The input variables ex ey ep and D are defined in planre The vector ed contains the nodal displacements a of the element and is obtained by the function extract as ed a u uz ug The output variables es o7 Oxz Oyy oz Ory ox oy et e Ezz Eyy Ezz Yay zz yz contain the stress and strain components The size of es and et follows the size of D Note that for plane stress 0 and for plane strain o 0 Theory The strains and stresses are computed according to e B a o De where the matrices D B and a are described in planre and where the evaluation point x y is chosen to be at the center of the element 5 9 15 ELEMENT plantce Two dimensional solid elements Purpose Compute element matrices for a rectangular Turner Clough element in plane strain or plane stress uz Us Xp 4 X33 A Ug U3 gt X Y2 Syntax Ke plantce ex ey ep Ke fe plantce ex ey ep eq Description plantce provides an element stiffness matrix Ke and an element load vector fe for a rectangular Turner Clough element in pla
22. plot x y plots vector x versus vector y Straight lines are drawn between each pair of values Various line types plot symbols and colors may be obtained with plot x y s where s is a 1 2 or 3 character string made from the following characters solid line point y yellow dotted line o circle m magenta dashdot line x x mark c cyan dashed line plus r red star g green b blue w white k black Default is solid blue lines Example The statement plot x y Xx y ro plots the data twice giving a solid blue line with red circles at the data points Note This is a MATLAB built in function For more information about the plot function type help plot 8 15 GRAPHICS print Purpose Create hardcopy output of current figure window Syntax print filename Description print with no arguments sends the contents of the current figure window to the default printer print filename creates a PostScript file of the current figure window and writes it to the specified file Note This is a MATLAB built in function For more information about the print function type help print GRAPHICS 8 16 text Purpose Add text to current plot Syntax text x y string Description text adds the text in the quotes to location x y on the current axes where x y is in units from the current plot If x and y are vectors text writes the text at all locations given If string
23. the spring stiffness in the axial direction k and the spring stiffness in the transverse direction ky The element load vector fe can also be computed if uniformly distributed loads are applied to the element The optional input variable eq described in beam2e contains the distributed loads 5 6 11 ELEMENT beam2w Two dimensional beam element Theory The element stiffness matrix K stored in Ke is computed according to Ke G K K G where G and K are described in beam2e The matrix K is given by 140k 0 0 70ka 0 0 156k kL 0 54k C 0 22kL 4k 0 13k L S 490 70ka 0 0 140k 0 0 54k 13kL 0 156k 0 13k 0 3k 0 22kL where the length L z 21 Y2 i ELEMENT 5 6 12 0 13k L 3k L 22k L Ak L Two dimensional beam element beam2ws Purpose Compute section forces in a two dimensional beam element on elastic foundation v2 X M lt N J Syntax es beam2ws ex ey ep ed es beam2ws ex ey ep ed eq Description beam2ws computes the section forces at the ends of the beam element on elastic foundation beam2w The input variables ex ey ep and eq are defined in beam2w and the element dis placements stored in ed are obtained by the function extract If distributed loads are applied to the element the variable eq must be included The output variable s Mn m T M V M contains the section forces at the ends of the beam 5 6
24. Assuming the matrix E was a zero matrix before the statement is executed the result will be I coc o woo oroo oo0oo0oo0o Note This is a MATLAB built in character MATRIX 3 4 Purpose Matrix arithmetic Syntax A B A B AxB A s Description Matrix operations are defined by the rules of linear algebra Examples An example of a sequence of matrix to matrix operations is D A B C A matrix to vector multiplication followed by a vector to vector subtraction may be defined by the statement b c Axx and finally to scale a matrix by a scalar s we may use B A s Note These are MATLAB built in operators 3 5 MATRIX abs Purpose Absolute value Syntax B abs A Description B abs A computes the absolute values of the elements of matrix A and stores them in matrix B Examples Assume the matrix 274 3 S The statement D abs C results in a matrix 7 4 ela stored in the workspace Note This is a MATLAB built in function For more information about the abs function type help abs MATRIX 3 6 det Purpose Matrix determinant Syntax a det A Description a det A computes the determinant of the matrix A and stores it in the scalar a Note This is a MATLAB built in function For more information about the det function type help det 3 7 MATRIX diag Purpose Diagonal matrices and diagonals of a matrix Syntax M diag v v diag
25. B VN where V 0 Oy Os 2 Oy Ox and where a a ar Ou Ox _ qTy 1 o lt E On ERI 5 ay Oy On Of On 5 0 21 ELEMENT plani4e Two dimensional solid elements If a larger D matrix than 3 x 3 is used for plane stress ptype 1 the D matrix is reduced to a 3 x 3 matrix by static condensation using oz Osz Oyz 0 These stress components are connected with the rows 3 5 and 6 in the D matrix respectively If a larger D matrix than 3 x 3 is used for plane strain ptype 2 the D matrix is reduced to a 3 x 3 matrix using Vez Yiz 0 This implies that a 3 x 3 D matrix is created by the rows and the columns 1 2 and 4 from the original D matrix Evaluation of the integrals is done by Gauss integration ELEMENT 5 5 22 Two dimensional solid elements plani4s Purpose Compute stresses and strains in a 4 node isoparametric element in plane strain or plane stress Us Syntax es et eci plani4s ex ey ep D ed Description plani4s computes stresses es and the strains et in a 4 node isoparametric element in plane strain or plane stress The input variables ex ey ep and the matrix D are defined in plani4e The vector ed contains the nodal displacements a of the element and is obtained by the function extract as ed a u U2 Us The output variables 1 1 1 1 1 1 Org Oyy oz Oxy ozz
26. BAd Coe teas ow Seon 5 2 1 5 3 Barmelements 8 22 a e a a ee Se ee ee Se A a 5 3 1 5 4 Heat flow elements 0 00 00 000 a a a a 5 4 1 5 9 Solid elements ese bed ish i ead eh He e aa E a a Se ee eed a 5 5 1 5 6 Beam elements a a a a a 5 6 1 Dep Plate element ea egal he er AWG arahna ena a Ue Pal betes ie an es 5 7 1 6 System functions 6 1 1 6 1 introductions ani ee Se a ees ee eke e Gide BO a eh oe i 6 1 1 6 2 Static systeni funtion s arp aa a a pe ote ee a alae d 6 2 1 6 3 Dynamic system functions ooo ae he gerd Gv ke ee Grn OS 6 3 1 7 Statements and macros 7 1 8 Graphics functions 8 1 9 User s Manual examples 9 1 1 9 troduction s s Ann aA AE A ak BOE A ek a Ee a ed 9 1 1 9 2 MATLAB introduction aoa aa a a 2 0 0 a a 9 2 1 Od Static nalysis sos sa be we oe a a Se ea ed ee ee k 9 3 1 Gels Dyn mic analysis s ei oh as aa Be Ges es a eee Dee ee De ee Te eT ae 9 4 1 9 5 Nonlinear analysis A ooa aa le hela A hed eed 9 5 1 1 Introduction The computer program CALFEM is a MATLAB toolbox for finite element applications This manual concerns mainly the finite element functions but it also contains descriptions of some often used MATLAB functions The finite element analysis can be carried out either interactively or in a batch oriented fashion In the interactive mode the functions are evaluated one by one in the MATLAB command window In the batch oriented mode a sequence of functions a
27. CALFEM diskette under the directory examples The example files are named according to the table 9 3 1 EXAMPLES exsl Static analysis Purpose Show the basic steps in a finite element calculation Description The general procedure in linear finite element calculations is carried out for a simple structure The steps are e define the model e generate element matrices e assemble element matrices into the global system of equations e solve the global system of equations e evaluate element forces Consider the system of three linear elastic springs and the corresponding finite element model The system of springs is fixed in its ends and loaded by a single load F 2k t VV AANA Vv Y The computation is initialized by defining the topology matrix Edof containing ele ment numbers and global element degrees of freedom gt gt gt Edof 1 1 2 2 2 3 3 23 the global stiffness matrix K 3x3 of zeros gt gt K zeros 3 3 K 0 0 0 0 0 0 0 0 0 EXAMPLES 9 3 2 Static analysis exsl and the load vector f 3x1 with the load F 100 in position 2 gt gt f zeros 3 1 2 100 Element stiffness matrices are generated by the function springle The element prop erty ep for the springs contains the spring stiffnesses k and 2k respectively where k 1500 gt gt k 1500 epl k ep2 2 k gt gt Kel springte ep1 Kel 1500 1500
28. Examples By the statement a 2 the scalar a is assigned a value of 2 An element in a matrix may be assigned a value according to A 2 5 3 The statement D 1 2 3 4 results in matrix 34 stored in the matrix bank To copy the contents of the matrix D to a matrix E use E D The character is used in the following statement to store the transpose of the matrix A in anew matrix F F A Note These are MATLAB built in characters MATRIX 3 2 Purpose Create vectors and do matrix subscripting Description The colon operator uses the following rules to create regularly spaced vectors j k isthe same as j j 1 k j i k is the same as j j i j 2i k The colon notation may also be used to pick out selected rows columns and elements of vectors and matrices A j is the j th column of A A i is the i th row of A Examples The colon used with integers d 1 4 results in a row vector d 1 2 3 4 stored in the workspace The colon notation may be used to display selected rows and columns of a matrix on the terminal For example if we have created a 3 times 4 matrix D by the statement D d 2xd 3 d resulting in 123 4 D 2 46 8 3 6 9 12 columns three and four are displayed by entering D 3 4 resulting in 3 4 D 3 4 6 8 9 12 In order to copy parts of the D matrix into another matrix the colon notation is used as E 3 4 2 3 D 1 2 3 4 3 3 MATRIX
29. M Description For a vector v with n components the statement M diag v results in an n x n matrix M with the elements of v as the main diagonal For an x n matrix M the statement v diag M results in a column vector v with n components formed by the main diagonal in M Note This is a MATLAB built in function For more information about the diag function type help diag MATRIX 3 8 full Purpose Convert sparse matrices to full storage class Syntax A full S Description A full S converts the storage of a matrix from sparse to full If A is already full full A returns A Note This is a MATLAB built in function For more information about the full function type help full 3 9 MATRIX inv Purpose Matrix inverse Syntax B inv A Description B inv A computes the inverse of the square matrix A and stores the result in the matrix B Note This is a MATLAB built in function For more information about the inv function type help inv MATRIX 3 10 length Purpose Vector length Syntax n length x Description n length x returns the dimension of the vector x Note This isa MATLAB built in function For more information about the length function type help length 3 11 MATRIX max Purpose Maximum element s of a matrix Syntax b max A Description For a vector a the statement b max a assigns the scalar b the maximum element of the vector
30. Nye 0 0 5 3 3 ELEMENT bar2g Two dimensional bar element The transformation matrix G contains the direction cosines amp 2 Ly _ Y2 y Nee Nyg L Nye Nag L where the length L V 22 1 ye m1 ELEMENT 5 3 4 Two dimensional bar element bar2s Purpose Compute normal force in a two dimensional bar element Syntax es bar2s ex ey ep ed Description bar2s computes the normal force in the two dimensional bar elements bar2e and bar2g The input variables ex ey and ep are defined in bar2e and the element nodal dis placements stored in ed are obtained by the function extract The output variable es N contains the normal force N Theory The normal force N is computed from EA N T 1 1 1 G af where E A L and the transformation matrix G are defined in bar2e The nodal displacements in global coordinates e T af u U2 U3 u are also shown in bar2e Note that the transpose of a is stored in ed 5 3 5 ELEMENT bar3e Three dimensional bar element Purpose Compute element stiffness matrix for a three dimensional bar element Us X 2 Z2 U3 Xp pZ Syntax Ke bar3e ex ey ez ep Description bar3e provides the element stiffness matrix Ke for a three dimensional bar element The input variables ex z x E A ep E A supply the element nodal coordinates 71 Y1 21 2 etc the modulus of elastici
31. T yz 2 2 2 2 2 2 es ao Ore Oyy oz Oxy or Oyz n2 n2 n2 n2 na n2 Orr Oyy lon Oxy lon Oyz 1 1 1 1 n1 Ere Eyy Ezz Yay Ge Ty T Y 2 2 2 2 2 2 E E E s T2 Y2 et eT Lx yy za Vey zz Vyz eci n2 n2 n2 n2 n2 n2 Tn2 Un2 Err Eyy Es Yay ee Vyz contain the stress and strain components and the coordinates of the integration points The index n denotes the number of integration points used within the ele ment cf plani4e The number of columns in es and et follows the size of D Note that for plane stress 0 and for plane strain o 0 5 9 23 ELEMENT plani4s Two dimensional solid elements Theory The strains and stresses are computed according to e Ba o De where the matrices D B and af are described in plani4e and where the integration points are chosen as evaluation points ELEMENT 5 5 24 Two dimensional solid elements plani4f Purpose Compute internal element force vector in a 4 node isoparametric element in plane strain or plane stress Syntax ef plani4f ex ey ep es Description plani4f computes the internal element forces ef in a 4 node isoparametric element in plane strain or plane stress The input variables ex ey and ep are defined in plani4e and the input variable es is defined in planiA4s The output variable ef f fa fiz x bd fis contains the components of the internal force vector Theory The internal force vector
32. The group of material functions comprises functions for constitutive models The available models can treat linear elastic and isotropic hardening von Mises material These material models are defined by the functions Material property functions hooke mises dmises Form linear elastic constitutive matrix Compute stresses and plastic strains for isotropic hardening von Mises material Form elasto plastic continuum matrix for isotropic hardening von Mises material MATERIAL hooke Purpose Compute material matrix for a linear elastic and isotropic material Syntax D hooke ptype E v Description hooke computes the material matrix D for a linear elastic and isotropic material The variable ptype is used to define the type of analysis plane stress plane strain axisymmetry three dimensional analysis ptype Hy WwW Dd The material parameters E and v define the modulus of elasticity F and the Poisson s ratio v respectively For plane stress ptype 1 D is formed as I v 0 E E v 1 0 ae ae l v For plane strain ptype 2 and axisymmetry ptype 3 D is formed as l v v V 0 E v l v v 0 i v 1 2 v ov i v 0 0 0 0 1 2 For the three dimensional case ptype 4 D is formed as l v v v 0 0 0 v l v vw 0 0 0 E v v l v 0 0 0 TAa 0 oO 0 10 2 0 0 0 0 0 0 1 2r 0 0 0 0 0 0 1 2v MATERIAL 4 2 mises Purpose Compute stresses and plastic st
33. a For a matrix A the statement b max A returns a row vector b containing the maximum elements found in each column vector in A The maximum element found in a matrix may thus be determined by c max max A Examples Assume the matrix B is defined as a oa aes The statement d max B results in a row vector d 3 4 The maximum element in the matrix B may be found by e max d which results in the scalar e 4 Note This is a MATLAB built in function For more information about the max function type help max MATRIX 3 12 min Purpose Minimum element s of a matrix Syntax b min A Description For a vector a the statement b min a assigns the scalar b the minimum element of the vector a For a matrix A the statement b min A returns a row vector b containing the min imum elements found in each column vector in A The minimum element found in a matrix may thus be determined by c min min A Examples Assume the matrix B is defined as 7 4 ae The statement d min B results in a row vector d 7 8 The minimum element in the matrix B is then found by e min d which results in the scalar e 8 Note This is a MATLAB built in function For more information about the min function type help min 3 13 MATRIX ones Purpose Generate a matrix of all ones Syntax A ones m n Description A ones m n results in an m times n matrix A with all ones Note This is a M
34. and the element load vector fe for a quadrilateral heat flow element The element nodal coordinates 71 Y1 2 etc are supplied to the function by ex and ey the element thickness t is supplied by ep and the thermal conductivities or corresponding quantities kre key etc are supplied by D Kiss ks ep t ea ex z T2 3 T4 ey y Y2 Y3 Ya yx yy If the scalar variable eq is given in the function the element load vector fe is com puted using eq Q where Q is the heat supply per unit volume Theory In computing the element matrices a fifth degree of freedom is introduced The location of this extra degree of freedom is defined by the mean value of the coordinates in the corner points Four sets of element matrices are calculated using flw2te These matrices are then assembled and the fifth degree of freedom is eliminated by static condensation ELEMENT 5 4 6 Two dimensional heat flow elements flw2qs Purpose Compute heat flux and temperature gradients in a quadrilateral heat flow element Syntax es et flw2qs ex ey ep D ed es et flw2qs ex ey ep D ed eq Description flw2qs computes the heat flux vector es and the temperature gradient et or corre sponding quantities in a quadrilateral heat flow element The input variables ex ey eq and the matrix D are defined in flw2qe The vector ed contains the nodal temperatures af of the element and is obtained by the function extract
35. applicable in modelling differ ent physical applications Table 3 below shows the relation between the primary variable a the constitutive matrix D and the load vector f for a chosen set of two dimensional physical problems Problem type a D f Designation Heat flow T Nag dg Q T temperature Az Ay thermal conductivity Q heat supply Groundwater flow kg ky Q piezometric head kz ky perme abilities Q fluid supply 1 1 St Venant torsion Go G zy Zu 20 stress function Gzy Gzz shear moduli angle of torsion per unit length Table 3 Problem dependent parameters 5 4 1 ELEMENT Heat flow elements T T T flw2qe flw2i4e flw2te flw2i8e flw3i8e 2D heat flow functions flw2te Compute element matrices for a triangular element flw2ts Compute temperature gradients and flux flw2qe Compute element matrices for a quadrilateral element flw2qs Compute temperature gradients and flux flw2i4e Compute element matrices 4 node isoparametric element flw2i4s Compute temperature gradients and flux flw2i8e Compute element matrices 8 node isoparametric element flw2i8s Compute temperature gradients and flux 3D heat flow functions flw3i8e Compute element matrices 8 node isoparametric element flw3i8s Compute temperature gradients and flux ELEMENT 5 4 2 Two dimensional heat flow eleme
36. diary function type help diary 2 3 GENERAL PURPOSE disp Purpose Display a variable in matrix bank on display screen Syntax disp A Description disp A displays the matrix A on the display screen Note This is a MATLAB built in function For more information about the disp function type help disp GENERAL PURPOSE 2 4 echo Purpose Control output on the display screen Syntax echo on echo off echo Description echo on turns on echoing of commands inside Script files echo off turns off echoing echo by itself toggles the echo state Note This is a MATLAB built in function For more information about the echo function type help echo 25 GENERAL PURPOSE format Purpose Control the output display format Syntax See the listing below Description format controls the output format By default MATLAB displays numbers in a short format with five decimal digits Command Result Example format short 5 digit scaled fixed point 3 1416 format long 15 digit scaled fixed point 3 14159265358979 format short e 5 digit floating point 3 1416e 00 format longe 16 digit floating point 3 141592653589793e 00 Note This is a MATLAB built in function For more information about the format func tion type help format GENERAL PURPOSE 2 6 help Purpose Display a description of purpose and syntax for a specific function Syntax help function name Description
37. ey ep Ed 3 EXAMPLES 9 3 10 Static analysis exs3 Results a Q 0 1 0e 03 0 0 0095 0 0 6 6667 0 0228 0 0000 0 0038 0 0 0 0000 0 0199 0 0000 0 0047 0 0 0 0000 0 0 0000 0 0076 0 3 3333 0 0000 esi 1 0e 04 0 0 6667 0 0000 0 0 6667 2 0000 es2 1 0e 04 x 0 0 3333 2 0000 0 0 3333 1 0000 es3 1 0e 04 x 0 0 3333 1 0000 0 0 3333 0 0000 9 3 11 EXAMPLES exs4 Static analysis Purpose Analysis of a plane truss Description Consider a plane truss loaded by a single force P 0 5 MN A 25 0e m E 2 10e MPa 2m 2m 30 The corresponding finite element model consists of ten elements and twelve degrees of freedom The topology is defined by the matrix gt gt Edof 1 o O aon AN ps ODNATBRWHN TOrRNWENN Owe eS N DMN CO FW FPN Onone ON Ww e Oo eS EXAMPLES 9 3 12 Static analysis exs4 A global stiffness matrix K and a load vector f are defined The load P is divided into x and y components and inserted in the load vector f gt gt K zeros 12 gt gt f zeros 12 1 11 0 5e6 sin pi 6 12 0 5e6 cos pi 6 The element matrices Ke are computed by the function bar2e These matrices are then assembled in the global stiffness matrix using the function assem gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt
38. help provides an online documentation for the specified function Example Typing gt gt help barie yields Ke bar1e ep PURPOSE Compute element stiffness matrix for spring analog element INPUT ep k spring stiffness or analog quantity OUTPUT Ke stiffness matrix dim Ke 2 x 2 Note This is a MATLAB built in function For more information about the help function type help help 2 GENERAL PURPOSE load Purpose Retrieve variable from disk and load in workspace Syntax load filename load filename eat Description load filename retrieves the variables from the binary file filename mat load filename ext reads the ASCII file filename ext with numeric data arranged in m rows and n columns The result is an m by n matrix residing in workspace with the name filename i e with the extension stripped Note This is a MATLAB built in function For more information about the load function type help load GENERAL PURPOSE 2 8 quit Purpose Terminate CALFEM session Syntax quit Description quit filename terminates the CALFEM without saving the workspace Note This is a MATLAB built in function For more information about the quit function type help quit 2 9 GENERAL PURPOSE Save Purpose Save workspace variables on disk Syntax save filename save filename variables save filename variables ascii Description save filename writes all variables resid
39. matrix Ed defined in extract If isov is a scalar it determines the number of iso lines to be displayed If isov is a vector it determines the values of the iso lines to be displayed number of iso lines equal to length of vector isov isov isolines isov isovalue 1 isovalue n The variable plotpar sets plot parameters for the iso lines plotpar linetype linecolor textfcn arrowtype 1 solid arrowcolor 1 black 2 dashed 2 blue 3 dotted 3 magenta 4 red textfcn 0 the iso values of the lines will not be printed 1 the iso values of the lines will be printed at the iso lines 2 the iso values of the lines will be printed where the cursor indicates Default is solid black lines and no iso values printed Limitations Supported elements are triangular 3 node and quadrilateral 4 node elements 8 9 GRAPHICS elprinc2 Purpose Draw element principal stresses as arrows for two dimensional elements Syntax magnfac elprinc2 Ex Ey Es magnfac elprinc2 Ex Ey Es plotpar elprinc2 Ex Ey Es plotpar magnfac Description elprinc2 displays element principal stresses for a number of elements of the same type The principal stresses are displayed as arrows at the element centroids Note that only the principal stresses are displayed To display the element mesh use eldraw2 Input variables are the coordinate matrices Ex and Ey and the element stresses matrix Es defined in plants or planqs Th
40. the vertex colors If X Y and C are matrices of the same size fill X Y C draws one polygon per column with interpolated colors Example The solution of a heat conduction problem results in a vector d with nodal tem peratures The temperature distribution in a group of triangular 3 node nen 3 or quadrilateral 4 node nen 4 elements with topology defined by edof can be displayed by Lex ey coordxtr edof Coord Dof nen ed extract edof d colormap hot fill ex ey ed Note This is a MATLAB built in function For more information about the fill function type help fill GRAPHICS 8 12 grid Purpose Grid lines for 2D and 3D plots Syntax grid on grid off grid Description grid on adds grid lines on the current axes grid off takes them off grid by itself toggles the grid state Note This is a MATLAB built in function For more information about the grid function type help grid 8 13 GRAPHICS hold Purpose Hold the current graph Syntax hold on hold off hold Description hold on holds the current graph hold off returns to the default mode where plot functions erase previous plots hold by itself toggles the hold state Note This is a MATLAB built in function For more information about the hold function type help hold GRAPHICS 8 14 plot Purpose Linear two dimensional plot Syntax plot x y plot x y linetype Description
41. variables are EXAMPLES 9 2 2 MATLAB introduction exil gt gt whos Name Size Bytes Class A 2x2 32 double array B 2x2 32 double array u 1x4 32 double array v 1x6 48 double array x 1x1 8 double array Grand total is 19 elements using 152 bytes The value of a variable is displayed by writing the variable name 1 2 3 4 and the dimension mxn of a variable is obtained by gt gt size u ans 1 4 where the answer is temporarily stored in the vector ans The variable x is removed from workspace by gt gt clear x To remove all variables from workspace clear without argument is used Assignment of a value to an element in a matrix is made as gt gt A 2 2 9 A 9 2 3 EXAMPLES exil MATLAB introduction To select a complete row or column colon notation is used gt gt s A 1 s 1 3 gt gt t A 2 t 3 9 A zero matrix K 4x4 is generated by gt gt K zeros 4 4 K oOo OO Oo oOoOO oOoOO 0 oO OO 0 Similarly an mxn matrix of all ones can be generated by ones m n Expand an already defined matrix gt gt H A B H Nowe on oO N gt gt J A B EXAMPLES 9 2 4 MATLAB introduction exi2 Purpose Show some examples of basic matrix and element by element operations Description Consider the following matrices a 5 12 3 b 10 4 1 6 3 254 ree Be oa The transpose of a matrix is generated by gt gt A ans 1 2 8 3 4
42. 0 1 0 1 1 nhistr 1 1 nev 4 tine integration parameters ip dt T 0 25 0 5 10 nev ntimes nhistr time integration Dsnapr Dr Vr Ar step2 kr mr dr0 vr0 ip fr mapping back to original coordinate system 9 4 9 EXAMPLES exd3 Dynamic analysis DsnapR Egv 1 nev Dsnapr DR Egv nhist 1 nev Dr 4 plot time history for two DOF s figure 1 plot t DR 1 t DR 2 axis 0 1 0000 0 0100 0 0200 grid xlabel time sec ylabel displacement m title Displacement time at the 4th and 11th gt degree of freedom text 0 3 0 017 solid line impact point x direction text 0 3 0 012 dashed line center horizontal beam y direction text 0 3 0 007 2 EIGENVECTORS ARE USED EXAMPLES 9 4 10 Dynamic analysis exd4 Purpose This example deals with a time varying boundary condition and time integration for the frame structure defined in exd1 Description Suppose that the support of the vertical beam is moving in the horizontal direction The commands below prepare the model for time integration Note that the structure of the boundary condition matrix bc differs from the structure of the load matrix f defined in exd2 The first eight eigenmodes Hz Time dependent boundary condition at the support DOF 1
43. 010 1 0005 0 0005 0 0014 1 0003 0 0005 0 0010 1 0005 0 0005 0 0014 es3 es3 499 5747 326 9323 981 6016 1 0e 03 499 5747 326 9323 979 9922 0 4996 0 4595 1 3793 0 4996 0 4595 1 3777 Using the second order theory the horizontal displacement of the upper left corner of the frame increases from 0 5 to 0 8 mm Both moments M4 and Mp are increased from 1 0 to 1 4 kNm 9 5 5 EXAMPLES exn2 Nonlinear analysis Purpose Buckling analysis of a plane frame Description The frame of exn1 is in this example analysed with respect to security against buckling for a case when all loads are increased proportionally The initial load distribution of exn1 is increased by a loading factor alpha until buckling occurs i e the deter minant of the stiffness matrix K passes zero For each value of alpha a second order theory calculation of type exn1 is performed The horizontal displacement a4 and the moment M4 are plotted against alpha The shape of the buckling mode is also plotted using the last computed displacement vector before buckling occurs The finite element model i e the vectors Edof ep1 ep3 Ex and Ey defined in exn1 is used 0 zeros 12 1 0 4 1000 0 5 1000000 0 8 1000000 j 0 for alpha 1 0 1 20 j jtt N 0 1 0 0 No 1 1 1 eps 0 00001 n 0 while abs N 1 NO 1 NO 1 gt eps n n 1 K zeros 12 12 f f0 alpha Kel beam2g Ex 1 Ey 1 ep1 N 1 Ke2 beam2g Ex 2
44. 12EI __6FI L 1 4 L 1 p L3 1 L 1 0 6EL 2EI 1 5 0 6BI 4BI 1 4 L2 1 p L p L2 1 p L p with LEI H TGAk ELEMENT 5 6 8 Two dimensional Timoshenko beam element beam2ts Purpose Compute section forces in a two dimensional Timoshenko beam element Syntax es beam2ts ex ey ep ed es beam2ts ex ey ep ed eq es edi eci beam2ts ex ey ep ed eq n Description beam2ts computes the section forces and displacements in local directions along the beam element beam2t The input variables ex ey ep and eq are defined in beam2t The element displace ments stored in ed are obtained by the function extract If distributed loads are applied to the element the variable eq must be included The number of evaluation points for section forces and displacements are determined by n If n is omitted only the ends of the beam are evaluated The output variables es NVM edi uv eci X consist of column matrices that contain the section forces the displacements and rotation of the cross section note that the rotation is not equal to 2 and the evaluation points on the local x axis The explicit matrix expressions are N V M i w A No Vo M U2 Uz b es f edi eci Nn Va Mn iin Um O pa where L is the length of the beam element Theory 5 6 9 ELEMENT beam2ts Two dimensional Timoshenko beam element The evaluation of the section forces is based on the solut
45. 2 magnfac Time num2str ntimes i text 3 i 1 5 1 5 Time end Eyt Ey 4 for i 6 10 Ext Ex i 6 3 eldraw2 Ext Eyt 2 3 0 Edb extract Edof Dsnap i eldisp2 Ext Eyt Edb 1 2 2 magnfac Time num2str ntimes i text 3 i 6 5 2 5 Time end Displacement time at the 1st 4th and 11th degree of freedom 0 02 i solid line bottom vertical beam x direction 0 027 7 dashed line center vertical beam x direction 0 015 dashed dotted line center horizontal beam y direction 0 005 0 005 0 01 0 015 v7 0 02 l Time history at DOF 1 DOF 4 and DOF 11 EXAMPLES 9 4 12 Dynamic analysis exd4 Snapshots sec magnification 20 _ k i i 0 3 k 0 4 1 0 5 X X 0 8 i 0 9 1 x x Snapshots of the deformed geometry for every 0 1 sec 9 4 13 EXAMPLES 9 5 Nonlinear analysis This section illustrates some nonlinear finite element calculations Nonlinear analysis exN1 Second order theory analysis of a frame exN2 Buckling analysis of a frame Note The examples listed above are supplied as m files on the CALFEM diskette under the directory examples The example files are named according to the table 9 5 1 EXAMPLES exnl Nonlinear analysis Purpose Analysis of a plane frame using second order theory Description The frame of exs5 is analysed
46. 33 D34 D35 D36 D Da Dog D 3 or D Day Da2 Daz Daa Das Dag D3 D32 D33 ELEMENT 5 5 4 Two dimensional solid elements plante If uniformly distributed loads are applied to the element the element load vector fe is computed The input variable by ata containing loads per unit volume by and b is then given Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to Ke C71 J B D BtdA C A C N btdA A with the constitutive matrix D defined by D and the body force vector b defined by eq The evaluation of the integrals for the triangular element is based on a linear dis placement approximation u x y and is expressed in terms of the nodal variables u1 Ug Ug as u z y N a N C a where Ug Jizy000 u N O0lary 1 Ti Yi 0 0 0 ui 00 0 1am U2 a 1 T2 Y2 0 0 0 e U3 Serger at esas a 1 T3 Y3 0 0 0 U5 0 0 0 1 T3 Y3 ug The matrix B is obtained as o 0 Ox 2g 2 yea o B VNc where V 0 gt Oy Oe ue Oy Ox If a larger D matrix than 3 x 3 is used for plane stress ptype 1 the D matrix is reduced to a 3 x 3 matrix by static condensation using oz Osz Oyz 0 These stress components are connected with the rows 3 5 and 6 in the D matrix respectively 5 9 9 ELEMENT plante Two dimensional solid elements If a larger D matrix than 3 x 3 is used for plane
47. 4 iterations are performed Determinant lt 0 buckling load is passed The requested plots of the horizontal displacement the moment M4 and the shape of the buckling mode are generated by the following commands figure 1 clf plot deform loadfact deform loadfact axis 0 0 1 1 3 5 grid xlabel Horizontal displacement m ylabel alpha title Displacement alpha for the upper left corner figure 2 clf plot M loadfact M loadfact axis 0 1e5 1 3 5 grid xlabel Moment in A Nm ylabel alpha title Supporting moment M A alpha figure 3 clf axis off eldraw2 Ex Ey 2 3 0 Edi extract Edof bmode eldisp2 Ex Ey Ed1 1 1 1 2 title Shape of buckling mode EXAMPLES 9 5 8 exn2 Nonlinear analysis Displacement alpha for the upper left corner 3 5 T T T 3h yee SSeS eS SS Ss Benin E SRO E E x Z 2 57 i 4 g 4 2 4 4 4 15H 4 n l L L L L i L L L L 0 0 01 0 02 0 03 0 04 0 05 0 06 0 07 0 08 0 09 0 1 horizontal displacement m Supporting moment M A alpha 3 5 T T T T T 3H Soaps ep ee Teme mmn m eas a ae 4 x 2 55 4 f E t 4 art 4 n 4 1 57 4 a 1 ii L L L 1 2 3 4 5 6 7 8 9 10 Moment in A Nm x 10 Shape of buckling mode 9 5 9 EXAMPLES BIDS ct ne oeior be 3 6 assem 6 2 2 AXIS
48. 5 0 025 0 05 0 025 0 0 05 0 025 0 05 0 05 0 05 0 0 075 0 025 0 075 0 05 0 075 0 0 1 0 025 0 1 0 05 0 1 J Dof 1 2 3 4 5 6 T 8 9 10 11 12 13 14 15 M TRER Element properties topology and coordinates ep 1 D 1 0 0 1 Edof 1 1 2 5 4 2 2 3 6 5 3 4 5 8 7 4 5 6 9 8 5 7 8 11 10 6 8 9 12 11 7 10 11 14 13 8 11 12 15 14 Ex Ey coordxtr Edof Coord Dof 4 Generate FE mesh eldraw2 Ex Ey 1 3 0 Edof 1 pause clf Create and assemble element matrices for i 1 8 Ke flw2qe Ex i Ey i ep D K assem Edof i K Ke end bc 1 0 2 0 3 0 4 0 7 0 10 0 13 0 5 e 3 14 1e 3 15 1e 3 EXAMPLES 9 3 30 Static analysis exs7 a Q solveq K f bc Ed extract Edof a for i 1 8 Es i flw2qs Ex i Ey i ep D Ed i end 4 Draw flux vectors and contour lines eldraw2 Ex Ey elflux2 Ex Ey Es pause clf eldraw2 Ex Ey eliso2 Ex Ey Ed 5 0 17 0 09 we 0 08 0 06 gt 0 05 0 03 0 01 F i L L i i 0 02 0 0 02 0 04 0 06 0 08 x Flux vectors 9 3 31 EXAMPLES exs7 Static analysis 0 07f gt 0 05F 0 04f 0 03 0 02 0 0 02 0 04 0 06 0 08 x Contour lines Two comments concerning the contour lines In the upper left corner the contour lines should physically have met at the corner point However the drawing of the contour l
49. 813165155695703e 02 2 916830902640915e 02 EXAMPLES 9 2 10 MATLAB introduction exid Purpose How to make a command window session log file Description The diary and echo commands are useful for presentation purposes since the complete set of statements can be saved together with some selected results The command gt gt diary filename saves all subsequent terminal input and resulting output in the file named filename on the default device The file is closed using gt gt diary off Consider the script file test m rane Script file test m diary testlog echo on A 0 4 2 8 B 3 9 5 7 C A B 2537 echo off diary off Normally the statements in an m file do not display during execution The com mands echo on and echo off allow the statements to be viewed as they execute Execution of test m yields gt gt test A 0 4 2 8 B 3 9 5 7 C A B 2537 C 0 0079 0 0110 0 0181 0 0292 in the command window and on the file testlog as well 9 2 11 EXAMPLES exi6 MATLAB introduction Purpose How to display vectors and handle the graphics window Description The contents of a vector versus the vector index or a vector versus another vector can be displayed in the graphics window Consider the vectors x 1 25 y 5 22 16 The function gt gt plot y plots the contents of the vector y versus vector index and gt gt plot x y plots the contents of the vector y versus the
50. ATLAB built in function For more information about the ones function type help ones MATRIX 3 14 red Purpose Reduce the size of a square matrix by omitting rows and columns Syntax B red A b Description B red A b reduces the square matrix A to a smaller matrix B by omitting rows and columns of A The indices for rows and columns to be omitted are specified by the column vector b Examples Assume that the matrix A is defined as 1 2 3 4 5 6 7 8 oie 9 10 11 12 13 14 15 16 and b as i The statement B red A b results in the matrix le B 5 d 3 15 MATRIX size Purpose Matrix dimensions Syntax d size A m n size A Description d size A returns a vector with two integer components d m n from the matrix A with dimensions m times n m n size A returns the dimensions m and n of the m x n matrix A Note This is a MATLAB built in function For more information about the size function type help size MATRIX 3 16 sparse Purpose Create sparse matrices Syntax S sparse A S sparse m n Description S sparse A converts a full matrix to sparse form by extracting all nonzero matrix elements If S is already sparse sparse S returns S S sparse m n generates an m times n sparse zero matrix Note This is a MATLAB built in function For more information about the sparse function type help sparse 3 17 MATRIX Spy Purpose Visualize matr
51. AX oenn 3 12 MIM rphng hi ae 3 13 MISES ieia cea 4 3 ONES psa orraa 3 14 plani4e 5 5 20 plani4f 5 5 25 plani4s 5 5 23 plani8e 5 5 26 plani8f 5 5 31 plani8s 5 5 29 plange Noa 24 5 5 9 plangs 5 5 11 planre 5 5 12 planrs 5 5 15 plantce 5 5 16 plantes 5 5 19 plante 5 5 4 plantf 2 23 v2 5 5 8 plants 5 5 7 platre 5 7 2 DIGS crnini 57 5 Plotka 8 15 PATT 42 2 eas 8 16 function T 5 10100 Ree area 7 6 GUNG pees acre 2 9 Tedia iaaa ee tot a 3 15 ia ATA 6 3 8 SAVE Ws deoleahna ers 2 10 SIZE i aan oe 3 16 solige a n 5 5 32 soli8f 5 5 37 SONGS n a 5 0 30 solveq 6 2 9 sparse Caw ed spectra 6 3 9 springle 5 2 4 springls 5 2 5 SPY n skews 3 18 SAPE tarana 3 19 statcon 6 2 10 Stepl ennei 6 3 10 STEP oasia 3 6 3 12 SUM arnoia deus 3 20 SWEEP Gaus cet es 6 3 14 TEXT atresia 8 17 UL sated acre 8 18 Ypes irate 2 11 what a a 2 12 while oh soe Bit 7 4 who nananana 2 13 xlabel 8 19 ylabel 8 19 ZETOS oonan 3 21 zlabel 8 19
52. Compute heat flux and temperature gradients in an 8 node isoparametric heat flow element Syntax es et eci flw2i8s ex ey ep D ed Description flw2i8s computes the heat flux vector es and the temperature gradient et or corre sponding quantities in an 8 node isoparametric heat flow element The input variables ex ey ep and the matrix D are defined in flw2i8e The vector ed contains the nodal temperatures af of the element and is obtained by the function extract as ed a T To T3 rO Ts The output variables de a Beek q q oe ie ag T T Ox Oy T T t Y a x et VT x Oy ei n n p nea r Oy contain the heat flux the temperature gradient and the coordinates of the integra tion points The index n denotes the number of integration points used within the element cf flw2i8e Theory The temperature gradient and the heat flux are computed according to VT B a q DVT where the matrices D B and a are described in flw2i8e and where the integration points are chosen as evaluation points 5 4 13 ELEMENT flw3i8e Three dimensional heat flow elements Purpose Compute element stiffness matrix for an 8 node isoparametric element Syntax Ke flw3i8e ex ey ez ep D Ke fe flw3i8e ex ey ez ep D eq Description flw3i8e provides the element stiffness conductivity matrix Ke and the element load vector fe for an 8 node isopara
53. Des Des If different D matrices are used in the Gauss points these D matrices are stored in a global vector D For numbering of the Gauss points see eci in plani8s D D2 D D2 ELEMENT 5 5 26 Two dimensional solid elements plani8e If uniformly distributed loads are applied to the element the element load vector fe is computed The input variable by containing loads per unit volume by and b is then given Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT DB tdA A f NT bia with the constitutive matrix D defined by D and the body force vector b defined by eq The evaluation of the integrals for the isoparametric 8 node element is based on a displacement approximation u 7 expressed in a local coordinates system in terms of the nodal variables u1 u2 U16 as ulg n N af where ui v E 0 NE 0O NE 0 o Eai sak NE 0 NS an 0 A FO U16 The element shape functions are given by NE IQA N 0 a Ng F0 90 G n NE ZAE N F0 904 m G n N il EVA n Ng F0 O0 a e n N The matrix B is obtained as o an 0 ss ia o B VN where V Qu Oy i Or Oy Ox 5 5 27 ELEMENT plani8e Two dimensional solid elements and where a a de de Ox _ qTy 1 Og o On a 2 oy ay Oy On E On If a larger D matrix than 3 x 3 is used for plane s
54. ES exs5 Static analysis EXAMPLES 9 3 22 Static analysis exs6 Purpose Set up a frame consisting of both beams and bars and illustrate the calculations by use of graphics functions Description A frame consists of horizontal and vertical beams and is stabilized with diagonal bars The frame with its coordinates and loading is shown in the left figure and the finite element model in the right In the following the statements for analysing the frame are given as an m file The matrices of the global system i e the stiffness matrix K the load vector f the coordinate matrix Coord and the corresponding degrees of freedom matrix Dof are defined by K zeros 18 18 f zeros 18 1 4 13 1 Coord 0 oror 9 3 23 EXAMPLES exs6 Static analysis Dof 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 The material properties the topology and the element coordinates for the beams and bars respectively are defined by 4 Element properties topology and coordinates epi 1 1 1 Edof1 1 1 2 3 7 8 9 2 7 8 9 13 14 15 3 4 5 6 10 11 12 4 10 11 12 16 17 18 5 T 8 9 10 11 12 6 13 14 15 16 17 18 Ex1 Ey1 coordxtr Edof1 Coord Dof 2 ep2 1 1 Edof2 7 1 2 10 11 8 7 8 16 17 9 7 8 4 5 10 13 14 10 11 Ex2 Ey2 coordxtr Edof2 Coord Dof 2 To check the model the finite element mesh can be drawn eldraw2 E
55. K Ke2 gt gt K assem Edof 3 K Ke3 K assem Edof 4 K Ke4 gt gt K assem Edof 5 K Ke5 gt gt bc 1 17 6 20 gt gt T solveq K f bc 17 0000 16 0912 15 5567 16 8995 17 6632 20 0000 The temperature values in the node points are given in the vector T After solving the system of equations the heat flow through the wall is computed using extract and spring1s gt gt edi extract Edof 1 T gt gt ed2 extract Edof 2 T gt gt ed3 extract Edof 3 T gt gt ed4 extract Edof 4 T gt gt edd extract Edof 5 T 9 3 7 EXAMPLES exs2 Static analysis gt gt qi spring1s ep1 ed1 qi 12 9825 gt gt q2 spring1s ep2 ed2 q2 12 9825 gt gt q3 spring1s ep3 ed3 q3 12 9825 gt gt q4 spring1s ep4 ed4 q4 12 9825 gt gt q5 spring1s ep5 ed5 q5 12 9825 The heat flow through the wall equal for all elements is q 13 0 W m EXAMPLES 9 3 8 Static analysis exs3 Purpose Analysis of a simply supported beam Description Consider the simply supported beam loaded by a single load f 10000 N applied at a point 1 meter from the left support The corresponding finite element mesh is also shown The following data apply to the beam Young s modulus FE 2 10e1 N m Cross section area A 45 3e7 m Moment of inertia T 2510e78 m f 10kN 9m 2 3 8 11 4 A A A gis 1 2 3 T
56. Nonlinear analysis The MATLAB introduction examples explain some basic concepts and introduce a set of standard MATLAB functions usable in the finite element context The static linear examples illustrate finite element analysis of different structures loaded by stationary loads The dynamic linear examples illustrate some basic features in dynamics such as modal analysis and time stepping procedures The examples of nonlinear analysis cover subjects such as second order theory and buckling 9 1 1 EXAMPLES 9 2 MATLAB introduction The examples in this section illustrate basic MATLAB concepts such as handling of workspace variables and functions The examples are MATLAB introduction exil Handling matrices exi2 Matrix and array operations exi3 Create and handle m files exi4 Display formats exi5 Create a session log file exi6 Graphic display of vectors 9 2 1 EXAMPLES exil MATLAB introduction Purpose Show how to create and handle matrices in MATLAB Description The following commands create a scalar x two vectors u and v and two matrices A and B Lines starting with the MATLAB prompt gt gt are command lines while the other lines show the results from these commands gt gt x 7 gt gt v 0 0 4 2 0 0 4000 0 8000 1 2000 1 6000 2 0000 gt gt A 1 2 3 4 A 1 2 4 gt gt B 5 6 7 8 5 6 T 8 Both brief and detailed listing of variables is possible gt gt who Your
57. The positive directions of qz qj and qz follow the local beam coordinate system The distributed torque qo is positive if directed in the local z direction i e from local y to local z All the loads are per unit length Theory The element stiffness matrix K is computed according to K GTK G where in which k EA and k k 0 12EIz L3 0 0 0 6Elz L2 0 12EIz L8 0 0 0 6EIz L2 KE 0 0 0 0 0 k 0 0 0 0 0 z i am Ses SF SS Ea N E SS OOS SS ses 0 0 0 0 k 0 0 0 0 Sh 9 BER yay 0 aoe 0 0 0 0 k 0 0 0 0 6Ely 4EIz F o 5 i o th L L2 0 0 0 0 k O 6E Iz 12ETz NG L3 L2 0 k 0 0 0 0 oe 0 es Oe 0 0 0 Oo 70 Aig Se ae and where Neg 0 0 0 0 0 i 0 0 0 DO 0 Nez 0 0 0 0 0 O Tse Nyz N O0 0 O Ney Nyy Nng 0 0 O nez Nyz Naz 0 0 OS GO e ies ite 0 0 0 0 Nag Nyy 0 0 0 O Nez Nyz 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 6 25 SO SO Oo amp oO oS SOO OS Oo oO oS oS 3 8g ee amp 8 z N XR ELEMENT beam3e Three dimensional beam element The element length L is computed according to L V e2 x1 Y2 y1 22 21 In the transformation matrix G ngg specifies the cosine of the angle between the x axis and 7 axis and so on The element load vector ff stored in fe is computed according to fe GIF where ELEMENT 5 6 26 Thr
58. a are to be prescribed the row number and the corresponding values are given in the boundary condition matrix dof Uy do u TOR fa 2 do fnbc Unbe where the first column contains the row numbers and the second column the corre sponding prescribed values If Q is given in the function reaction forces are computed according to Q Ka f 6 2 9 SYSTEM statcon Static system functions Purpose Reduce system of equations by static condensation Syntax K1 f1 statcon K f b Description statcon reduces a system of equations Ka f by static condensation The degrees of freedom to be eliminated are supplied to the function by the vector dof dof dof nb where each row in b contains one degree of freedom to be eliminated The elimination gives the reduced system of equations Ki ay fi where K and f are stored in K1 and f1 respectively SYSTEM 6 2 10 6 3 Dynamic system functions The group of system functions comprises functions for solving linear dynamic systems by time stepping or modal analysis functions for frequency domain analysis etc Dynamic system functions dyna2 dyna2f fft freqresp gfunc ifft ritz spectra step1 step2 sweep Solve a set of uncoupled second order differential equations Solve a set of uncoupled second order differential equations in the frequency domain Fast Fourier transform Compute frequency response Linear interpolation between equally sp
59. aced points Inverse Fast Fourier transform Compute approximative eigenvalues and eigenvectors by the Lanc zos method Compute seismic response spectra Carry out step by step integration in first order systems Carry out step by step integration in second order systems Compute frequency response function Note Eigenvalue analysis is performed by using the function eigen see static system functions 6 3 1 SYSTEM dyna2 Dynamic system functions Purpose Compute the dynamic solution to a set of uncoupled second order differential equa tions Syntax X dyna2 w2 xi f g dt Description dyna2 computes the solution to the set i 2 amp witi wrx fig t STs ce of differential equations where g t is a piecewise linear time function The set of vectors w2 xi and f contains the squared circular frequencies w the damping ratios and the applied forces f respectively The vector g defines the load function in terms of straight line segments between equally spaced points in time This function may have been formed by the command gfunc The dynamic solution is computed at equal time increments defined by dt Including the initial zero vector as the first column vector the result is stored in the m by n matrix X n 1 being the number of time steps Note The accuracy of the solution is not a function of the output time increment dt since the command produces the exact solution for straight line
60. again but it is now subjected to a load case including a horizontal load and two vertical point loads Second order theory is used me i eres 5 A L E 4 0m A Ap I E Apl E Je Ay A B lh 6 0m Ag I E P H The finite element model consisting of three beam elements freedom is repeated here Edof 1 123456 2 10 11 12789 3 456789 EXAMPLES 9 5 2 45 3e m 2510e78 mt 142 8e74 m 33090e78 m4 2 10e gt MPa 1 MN 0 001 MN and twelve degrees of Nonlinear analysis exnl cae Element properties and global coordinates F 2 1e11 A1 45 3e 4 A2 142 8e 4 I1 2510e 8 I2 33090e 8 ep1 E A1 I1 ep3 E A2 12 Ex 0 0 6 6 0 6 Ey 0 4 0 4 4 4 f zeros 12 1 f 4 1000 5 1000000 8 1000000 The beam element function of the second order theory beam2g requires a normal force as input variable In the first iteration this normal force is chosen to zero This means that the first iteration is equivalent to a linear first order analysis using beam2e Since the normal forces are not known initially an iterative procedure has to be applied where the normal forces N are updated according to the results of the former iteration The iterations continue until the difference in normal force of the two last iteration steps is less than an accepted error eps N NO NO lt eps The small value given to the initial normal force N 1 is to avoi
61. and capacity matrices respectively The initial conditions are given by the vector d0 containing initial values of d The time integration procedure is governed by the parameters given in the vector ip defined as ip dt Ta nsnap nhist time dofi oo list of list of nsnap nhist moments dofs where dt specifies the time increment in the time stepping scheme T total time and a a time integration constant see 1 The parameter nsnap denotes the number of snapshots stored in Tsnap The selected elapsed times are specified in time whereas nhist is the number of time histories stored in D and V The selected degrees of freedom are specified in dof The following table lists frequently used values of a a 0 Forward difference forward Euler a 5 Trapezoidal rule Crank Nicholson a 1 Backward difference backward Euler The matrix f contains the time dependent load vectors If no external loads are active the matrix corresponding to f should be replaced by The matrix f contains the time dependent prescribed values of the field variable If no field variables are prescribed the matrix corresponding to pbound should be replaced by Matrix f is organized in the following manner SYSTEM 6 3 10 Dynamic system functions step1 time history of the load at dof time history of the load at dof time history of the load at dofn The dimension of f is number of degrees of freedom x
62. as ed a Ti To T3 T The output variables es q qs qy OT OT et VT Sa Ox Oy contain the components of the heat flux and the temperature gradient computed in the directions of the coordinate axis Theory By assembling four triangular elements as described in flw2te a system of equations containing 5 degrees of freedom is obtained From this system of equations the unknown temperature at the center of the element is computed Then according to the description in flw2ts the temperature gradient and the heat flux in each of the four triangular elements are produced Finally the temperature gradient and the heat flux of the quadrilateral element are computed as area weighted mean values from the values of the four triangular elements If heat is supplied to the element the element load vector eq is needed for the calculations 5 4 7 ELEMENT flw2i4e Two dimensional heat flow elements Purpose Compute element stiffness matrix for a 4 node isoparametric heat flow element T T x4Y4 3 93 Xpy X Y2 Syntax Ke flw2i4e ex ey ep D Ke fe flw2i4e ex ey ep D eq Description flw2i4e provides the element stiffness conductivity matrix Ke and the element load vector fe for a 4 node isoparametric heat flow element The element nodal coordinates x1 Y1 2 etc are supplied to the function by ex and ey The element thickness t and the number of Gauss points n
63. atrix is generated by as sembling this element stiffness matrix Ke in the global stiffness matrix K gt gt K assem Edof K Ke Finally the solution is calculated by defining the boundary conditions bc and solving the system of equations The boundary condition at dof 13 is set to 0 5x10 as an average of the concentrations at the neighbouring boundaries Displacements a and reaction forces Q are computed by the function solveq gt gt bc 1 0 2 0 3 0 4 0 7 0 10 0 13 0 5e 3 14 1e 3 15 1e 3 gt gt a Q solveq K f bc The element flows q are calculated from element concentration Ed gt gt Ed extract Edof a gt gt for i 1 8 gt gt es flw2qs ex ey ep D Ed i gt gt end EXAMPLES 9 3 28 Static analysis exs7 Results a Q 1 0e 03 1 0e 03 0 0 0165 0 0 0565 0 0 0399 0 0 0777 0 0662 0 0000 0 0935 0 0 0 2143 0 1786 0 0000 0 2500 0 0000 0 0 6366 0 4338 0 0000 0 5494 0 0000 0 5000 0 0165 1 0000 0 7707 1 0000 0 2542 Es 0 0013 0 0013 0 0005 0 0032 0 0049 0 0022 0 0020 0 0054 0 0122 0 0051 0 0037 0 0111 0 0187 0 0213 0 0023 0 0203 9 3 29 EXAMPLES exs7 Static analysis The following m file shows an alternative set of commands to perform the diffusion analysis of exs7 By use of global coordinates an FE mesh is generated Also plots with flux vectors and contour lines are created K zeros 15 f zeros 15 1 Coord 0 0 0 025 0 0 05 0 0 0 025 0 02
64. beam is 3 m and of the horizontal beam 2 m The following data apply to the beams vertical beam horizontal beam Young s modulus N m 3e10 3e10 Cross section area m 0 1030e 2 0 0764e 2 Moment of inertia m 0 171e 5 0 0801e 5 Density kg m 2500 2500 2m O S A a b The structure is divided into 4 elements The numbering of elements and degrees of freedom are apparent from the figure The following m file defines the finite element model material data E 3e10 rho 2500 Av 0 1030e 2 Iv 0 0171e 4 IPE100 Ah 0 0764e 2 Th 0 00801e 4 IPE80 epv E Av Iv rho Av eph E Ah Ih rhoxAh EXAMPLES 9 4 2 Dynamic analysis exdl h topology Edof 1 1 2 3 4 5 6 2 4 5 6 789 3 7 8 9 10 11 12 4 10 11 12 13 14 15 list of coordinates Coord 0 0 O 1 5 O 3 1 3 2 3 i list of degrees of freedom Dof 1 2 3 456 789 10 11 12 13 14 15 generate element matrices assemble in global matrices K zeros 15 M zeros 15 Ex Ey coordxtr Edof Coord Dof 2 for i 1 2 k m c beam2d Ex i Ey i epv K assem Edof i K k M assem Fdof i M m end for i 3 4 k m c beam2d Ex i Ey i eph K assem Edof i K k M assem Fdo
65. bles in workspace whos Working with files and controlling the command window diary Save session in a named file echo Control output on the display screen format Control the output display format quit Stop execution and exit from the CALFEM program 2 1 GENERAL PURPOSE clear Purpose Remove variables from workspace Syntax clear clear namel name names Description clear removes all variables from workspace clear namel name name8 removes specified variables from workspace Note This is a MATLAB built in function For more information about the clear function type help clear GENERAL PURPOSE 222 diary Purpose Save session in a disk file Syntax diary filename diary off diary on Description diary filename writes a copy of all subsequent keyboard input and most of the resulting output but not graphs on the named file If the file filename already exists the output is appended to the end of that file diary off stops storage of the output diary on turns it back on again using the current filename or default filename diary if none has yet been specified The diary function may be used to store the current session for later runs To make 2 this possible finish each command line with semicolon to avoid the storage of intermediate results on the named diary file Note This is a MATLAB built in function For more information about the
66. d division by zero in the second convergence check If N does not converge in 20 steps the analysis is interrupted Initial values for the iteration eps 0 01 h Error norm N 0 01 0 O Initial normal forces No 1 1 1 Normal forces of the initial former iteration n 0 Iteration counter eae Iteration procedure while abs N 1 NO 1 NO 1 gt eps n n 1 K zeros 12 12 Kel beam2g Ex 1 Ey 1 ep1 N 1 Ke2 beam2g Ex 2 Ey 2 ep1 N 2 Ke3 beam2g Ex 3 Ey 3 ep3 N 3 K assem Edof 1 K Kel K assem Edof 2 K Ke2 K assem Edof 3 K Ke3 9 5 3 EXAMPLES exnl Nonlinear analysis bc 1 0 2 0 3 0 10 0 11 0 12 O a solveq K f bc Ed extract Edof a esi beam2gs Ex 1 Ey 1 ep1 Ed 1 N 1 es2 beam2gs Ex 2 Ey 2 ep1 Ed 2 N 2 es3 beam2gs Ex 3 Ey 3 ep3 Ed 3 N 3 NO N N esi 1 1 es2 1 1 es3 1 1 if n gt 20 disp The solution doesn t converge return end end Displacements and element forces from the linear elastic analysis and from the second order theory analysis respectively EXAMPLES 0005 0042 0000 0005 0042 0000 9 5 4 0008 0042 0000 0008 0042 0000 Nonlinear analysis exnl esi esi 1 0e 05 1 0e 05 9 9967 0 0050 0 0102 9 9954 0 0050 0 0142 9 9967 0 0050 0 0098 9 9954 0 0050 0 0138 es2 es2 1 0e 06 1 0e 06 1 0003 0 0005 0 0
67. dV V with the constitutive matrix D defined by D and the body force vector b defined by eq The evaluation of the integrals for the isoparametric 8 node solid element is based on a displacement approximation u 7 expressed in a local coordinates system in terms of the nodal variables u1 us U24 as u f 7 N a where ul uz Ne 0 0 NE 0 0 NE 0 0 ie u u N 0 NE 0 0 NES 0 0 NE O a uz 0 0 NE 0 0 N 0 0 NE u24 The element shape functions are given by N 50 90 Md 9 NE 0 9A M0 9 n iieii iOa N Z04904 0 MESHE Ng 0 00 NE EA The B matrix is obtained as Be VN 5 9 33 ELEMENT soli8e Three dimensional solid elements where o T 0 0 o oy a a V i Se le o oO 0 Oy On Oy Ox o o Of 0 a Oz Ox o o ot On a ja a a dy E ot On O E On O Evaluation of the integrals is done by Gauss integration ELEMENT 5 5 34 Three dimensional solid elements soli8s Purpose Compute stresses and strains in an 8 node isoparametric solid element Syntax es et eci soli8s ex ey ez ep D ed Description soli8s computes stresses es and the strains et in an 8 node isoparametric solid element The input variables ex ey ez ep and the matrix D are defined in soli8e The vector ed contains the nodal displacements af of the element and is obtained by the function extract as ed a u U2 uz The output variables
68. default working directory is according to initial settings for example C USER A new working directory can be chosen by typing for example gt gt cd A which makes the root directory in drive A the working directory Files containing MATLAB and or CALFEM statements are characterized by the suffix m For example the file bar2e m contains statements to evaluate the element stiffness matrix for a two dimensional bar element An m file is allowed to include references to other m files and also to call itself recursively Two types of m files can be used script files and function files Script files collect a sequence of statements under one command name Function files allow new functions with input and or output variables to be defined Both script files and function files are ordinary ASCII text files and can therefore be created in an arbitrary editor In the MATLAB environment an m file editor can be activated from the pull down menu on top of the MATLAB window An example of a script file is given below The following sequence of statements is typed in the m file editor and saved as test m he ae Script file test m A 0 4 2 8 B 3 9 5 7 C A B V ea ea end A line starting with an is regarded as a comment line The statements are executed by writing the file name without the suffix m in the command window gt gt test A 0 4 2 8 EXAMPLES 9 2 8 MATLAB introduction exi3 B
69. described in beam2e and K is given by re 0 o 0 wes Th 0 gew 9 we Te 20 BA 0 0 EA L L e 0 0 eb by 0 For axial compression N lt 0 we have kL kL 1 get ieot a aE PRT eo P g 1 3 da z T z Ps 102 with T k rv For axial tension N gt 0 we have kL kL 1 kL coth Aaa ee a 2 2 ERA 1 3 da 571 3 Ps 102 with T k Vp The parameter p is given by _ NE P PEI ELEMENT 5 6 16 0 0 Hih Hio SElp g 0 0 121 SEL E ha te 3 3 To te o 3 To o Two dimensional beam element beam2g The equivalent nodal loads ff stored in the variable fe are computed according to fe GT fe where 1 T _ i i re Tk o S E a fo 5 TY 3 For an axial compressive force we have TN 2 _ l coskL L kLsinkL and for an axial tensile force eme 2 kLsinhkL kL 5 6 17 ELEMENT beam2gs Two dimensional beam element Purpose Compute section forces in a two dimensional nonlinear beam element v2 X M lt N J Syntax es beam2gs ex ey ep ed N es beam2gs ex ey ep ed N eq Description beam2gs computes the section forces at the ends of the nonlinear beam element beam2g The input variables ex ey and ep are defined in beam2e and the variables N and eq in beam2g The element displacements stored in ed are obtained by the function extract If a distributed transversal load is applied to the element the variable eq m
70. e loop sequence The different relation operators that can be used can be found under the if command Examples A loop continuing until a equals b while a b end Note This is MATLAB built in language STATEMENTS 7 4 function Purpose Define a new function Syntax function out1 out2 J name inl in2 Description name is replaced by the name of the function The input variables inl in2 can be scalars vectors or matrices and the same holds for the output variables out1 out2 Example To define the CALFEM function springle a file named springle m is created The file contains the following statements function Ke spring1le k Define the stiffness matrix for a one dimensional spring with spring stiffness k Ke k k k k i e the function springle is defined to return a stiffness matrix Note This is MATLAB built in language 7 5 STATEMENTS script Purpose Execute a stored sequence of statements Syntax name Description name is replaced by the name of the script Example The statements below are stored in a file named spring m and executed by typing spring in the MATLAB command window Stiffness matrix for a one dimensional spring with stiffness k 10 k 10 Ke springle k Note This is MATLAB built in language STATEMENTS T 6 8 Graphics functions The group of graphics functions comprises functions for elemen
71. e variable plotpar sets plot parameters for the principal stress arrows plotpar arrowtype arrowcolor arrowtype 1 solid arrowcolor 1 black 2 dashed 2 blue 3 dotted 3 magenta 4 red Default if plotpar is omitted is solid black arrows The magnification factor magnfac is a scalar that magnifies the arrows in relation to the element size The magnification factor is set automatically if it is omitted in the input list Limitations Supported elements are triangular 3 node and quadrilateral 4 node elements GRAPHICS 8 10 figure Purpose Create figures graph windows Syntax figure h Description figure h makes the h th figure the current figure for subsequent plotting functions If figure h does not exist a new one is created using the first available figure handle Note This isa MATLAB built in function For more information about the figure function type help figure 8 11 GRAPHICS fill Purpose Filled 2D polygons Syntax fill x y c fill X Y C Description fill x y c fills the 2D polygon defined by vectors x and y with the color specified by c The vertices of the polygon are specified by pairs of components of x and y If necessary the polygon is closed by connecting the last vertex to the first If c is a vector of the same length as x and y its elements are used to specify colors at the vertices The color within the polygon is obtained by bilinear interpolation in
72. ed in beam2e and the element displace ments stored in ed are obtained by the function extract If distributed loads are applied to the element the variable eq must be included The number of evaluation points for section forces and displacements are determined by n If n is omitted only the ends of the beam are evaluated The output variables es N VM edi ut 7V eci x consist of column matrices that contain the section forces the displacements and the evaluation points on the local x axis The explicit matrix expressions are Ni Vi M Uy V4 Nz V2 Me Uz Vz g es edi f eci f ak PR Na Va Ma Un Un L where L is the length of the beam element Theory 5 6 5 ELEMENT beam2s Two dimensional beam element The evaluation of the section forces is based on the solutions of the basic equations au dtu _4 9 0 ye age 4 de From these equations the displacements along the beam element are obtained as the sum of the homogeneous and the particular solutions a _ a a utu where gz Lz x lip Z BA 7 u N C Ga ty 3 dp ag he e Ti 24E1I L and 1000 0 0 0 010 0 8 U1 1z00 0 0 0001 0 0 re a err E s E 0O0O1LEP B Pe 0 0 0 1 2L 3r The transformation matrix G and nodal displacements af are described in beam2e Note that the transpose of a is stored in ed Finally the section forces are obtained from d d o gu N EA V EI M EI dz dz dx Exa
73. ee dimensional beam element beam3s Purpose Compute section forces in a three dimensional beam element xM y n Syntax es beam3s ex ey ez eo ep ed es beam3s ex ey ez eo ep ed eq es edi eci beam3s ex ey ez eo ep ed eq n Description beam3s computes the section forces and displacements in local directions along the beam element beam3e The input variables ex ey ez eo and ep are defined in beam3e and the element displacements stored in ed are obtained by the function extract If distributed loads are applied to the element the variable eq must be included The number of evaluation points for section forces and displacements are determined by n If n is omitted only the ends of the beam are evaluated The output variables es N V V T M M edi u w eci x consist of column matrices that contain the section forces the displacements and the evaluation points on the local z axis The explicit matrix expressions are Ni Va Vex T May W No Vp Va T Mp Mz es r 7 Nn Vin Ven T Myn Mzn 5 6 27 ELEMENT beam3s Three dimensional beam element 0 Uy Uy Wi YI T2 f Ug U2 W2 2 f edi i i eci ai Un Un Wn Pn L where L is the length of the beam element Theory The evaluation of the section forces is based on the solutions of the basic equations du dtu EA z 0 El 0 de 4 a Y d w d o EI as 0 GK Q 0 vag qe 4 F
74. ee node and quadrilateral four node elements 8 5 GRAPHICS eldia2 Purpose Draw the section force diagrams of a two dimensional beam element Ya wa TN Syntax magnfac eldia2 ex ey es eci eldia2 ex ey es eci magnfac eldia2 ex ey es eci magnfac magnitude Description eldia2 plots a section force diagram of a two dimensional beam element in its global position The input variables ex and ey are defined in beam2e and the input variables Sy T S2 T2 es eci f Sn n consist of column matrices that contain section forces and corresponding local x coordinates respectively The values in es and eci are computed in beam2s It should be noted however that whereas all three section forces are computed in beam2s only one of them shall be given as input to eldia2 by es GRAPHICS 8 6 eldia2 The magnification factor magnfac is a scalar that magnifies the diagram for visibility magnfac is set automatically if it is omitted in the input list The optional input magnitude S x y adds a scaled bar with length equivalent to a reference force S starting at coordi nates x y If no coordinates are given the starting point will be 0 0 5 Limitations Supported elements are two dimensional beam elements 8 7 GRAPHICS elflux2 Purpose Draw element flow arrows for two dimensional elements Syntax magnfac elflux2 Ex Ey Es magnfac elflux2 Ex Ey Es plotpar elflu
75. em functions coordxtr Purpose Extract element coordinates from a global coordinate matrix Ex x X3 X7 Xg Ey y Y2 Y7 Ye Syntax Ex Ey Ez coordxtr Edof Coord Dof nen Description coordxtr extracts element nodal coordinates from the global coordinate matrix Coord for elements with equal numbers of element nodes and dof s Input variables are the element topology matrix Edof defined in assem the global coordinate matrix Coord the global topology matrix Dof and the number of element nodes nen in each element T Y 2 k ty ea m 2 Y2 22 ko lb m Coord T3 Y3 23 Dof ks Ig ms nen nen The nodal coordinates defined in row 7 of Coord correspond to the degrees of freedom of row 2 in Dof The components k l and m define the degrees of freedom of node i and n is the number of global nodes for the considered part of the FE model 6 2 3 SYSTEM coordxtr Static system functions The output variables Ex Ey and Ez are matrices defined according to 1 1 1 1 Ly v2 T3 Ae Lnen 2 2 2 2 Ly 2 T3 Das XLnen Ex nel nel nel nel T T2 T3 e Lnen where row 7 gives the x coordinates of the element defined in row i of Edof and where nel is the number of considered elements The element coordinate data extracted by the function coordxtr can be used for plotting purposes and to create input data for the element stiffness functions SYSTEM 6 2 4 Static sys
76. ent variables and their size Examples who Your variables are A B C K M X k lambda whos name size elements bytes density complex A 3 by 3 9 72 Full No B 3 by 3 9 72 Full No C 3 by 3 9 72 Full No K 20 by 20 400 3200 Full No M 20 by 20 400 3200 Full No X 20 by 20 400 3200 Full No k 1 by 1 1 8 Full No lambda 20 by 1 20 160 Full No Grand total is 1248 elements using 9984 bytes Note These are MATLAB built in functions For more information about the functions type help who or help whos 2 13 GENERAL PURPOSE Purpose Continuation Syntax Description An expression can be continued on the next line by using Note This is a MATLAB built in function GENERAL PURPOSE 2 14 Purpose Write a comment line Syntax arbitrary text Description An arbitrary text can be written after the symbol Note This is a MATLAB built in character 2 15 GENERAL PURPOSE 3 Matrix functions The group of matrix functions comprises functions for vector and matrix operations and also functions for sparse matrix handling MATLAB has two storage modes full and sparse Only nonzero entries and their indices are stored for sparse matrices Sparse matrices are not created automatically But once initiated sparsity propagates Operations on sparse matrices produce sparse matrices and operations on a mixture of sparse and full matrices also normally produce sparse matrices The f
77. ep1 Ed 2 eq2 20 es2 1 0e 05 2 2533 0 1613 0 2185 2 2533 0 1613 0 1845 2 2533 0 1613 0 4266 gt gt es3 beam2s ex3 ey3 ep3 Ed 3 eq3 20 es3 1 0e 05 0 1613 2 2467 0 4070 0 1613 2 0099 0 2651 0 1613 2 2533 0 4266 EXAMPLES 9 3 18 Static analysis exs5 Section force diagrams are displayed using the function eldia2 gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt gt figure 1 magnfac eldia2 exl ey1 est 1 ecil1 magnitude 3e5 0 5 0 eldia2 ex1 ey1 esi 1 ecil magnfac eldia2 ex2 ey2 es2 1 eci2 magnfac eldia2 ex3 ey3 es3 1 eci3 magnfac magnitude axis 1 5 7 0 5 5 5 figure 2 magnfac eldia2 ex3 ey3 es3 2 eci3 magnitude 3e5 0 5 0 eldia2 ex1 ey1 es1 2 ecil magnfac eldia2 ex2 ey2 es2 2 eci2 magnfac eldia2 ex3 ey3 es3 2 eci3 magnfac magnitude axis 1 5 7 0 5 5 5 figure 3 magnfac eldia2 ex3 ey3 es3 3 eci3 magnitude 3e5 0 5 0 eldia2 ex1 ey1 es1 3 ecil magnfac eldia2 ex2 ey2 es2 3 eci2 magnfac eldia2 ex3 ey3 es3 3 eci3 magnfac magnitude axis 1 5 7 0 5 5 5 9 3 19 EXAMPLES exs5 Static analysis OF 300000 EXAMPLES 9 3 20 Static analysis exs5 OF k 300000 9 3 21 EXAMPL
78. esistance m 0 18 m2 K W concrete A 1 7 W mK thermal insulation A 0 04 W mK concrete A 1 7 W mK surface thermal resistance m 0 07 m2 K W 17 C Ty 0 070 m 0 100 m 0 100 m T T T T T Ts SO SSS SO 1 2 3 4 5 The wall is subdivided into five elements and the one dimensional spring analogy element springle is used Equivalent spring stiffnesses are k AA L for thermal conductivity and k A m for thermal surface resistance Corresponding spring stiffnesses per m of the wall are ky 1 0 07 14 2857 W C ka 1 7 0 070 24 2857 W C k3 0 040 0 100 0 4000 W C ka 1 7 0 100 17 0000 W C ks 1 0 18 5 5555 W C A global stiffness matrix K and a load vector f are defined The element matrices Ke are computed using springle and the function assem assembles the global stiffness matrix The system of equations is solved using solveq with considerations to the boundary conditions in bc The prescribed temperatures are T 17 C and T 20 C EXAMPLES 9 3 6 Static analysis exs2 gt gt Edof 1 2 oP WN EF Oo FW NY 3 4 5 gt gt K zeros 6 gt gt f zeros 6 1 gt gt ep1 1 0 07 ep2 1 7 0 07 gt gt ep3 0 040 0 10 ep4 1 7 0 10 gt gt ep5 1 0 18 gt gt Kel springte ep1 Ke2 springie ep2 gt gt Ke3 springle ep3 Ke4 springte ep4 gt gt Ke5 springle ep5 gt gt K assem Edof 1 K Ke1 K assem Edof 2
79. ever the response as a function of frequency can be computed efficiently In particular it is possible to identify the participating frequencies 6 3 5 SYSTEM gfunc Dynamic system functions Purpose Form vector with function values at equally spaced points by linear interpolation 4 e t tp8 t4 t 8 t2 t3 t3 P ty 8 ty tp8 t t5 2 ts Syntax t g gfunc G dt Description gfunc uses linear interpolation to compute values at equally spaced points for a discrete function g given by ty g t t g ta ty g tn as shown in the figure above Function values are computed in the range t lt t lt ty at equal increments dt being defined by the variable dt The number of linear segments steps is ty t1 dt The corresponding vector t is also computed The result can be plotted by using the command plot t g SYSTEM 6 3 6 Dynamic system functions ifft Purpose Transform function in frequency domain to time domain Syntax x ifft y x ifft y N Description ifft transforms a function in the frequency domain to a function in the time domain The function to be transformed is given in the vector y Each row in y contains Fourier terms in the interval oo lt w lt 00 The fft command can be used with one or two input arguments The scalar variable N can be used to specify the numbers of frequencies used in the Fourier transform The size of the o
80. expression else aid Description if initiates a conditional jump If logical expression produces the value True the statements following if are executed If logical expression produces the value False the next conditional statement elseif is checked elseif works like if One or more of the conditional statement elseif can be added after the initial conditional statement if If else is present the statements following else are executed if the logical expressions in all if and elseif statements produce the value False The if loop is closed by end to define the loop sequence The following relation operators can be used equal gt greater than or equal to gt greater than lt less than or equal to lt less than not equal Note This is MATLAB built in language STATEMENTS T 2 for Purpose Initiate a loop Syntax for i start inc stop end Description for initiates a loop which terminates when i gt stop The for loop is closed by end to define the loop sequence Examples fri 1 10 i takes values from 1 to 10 fri 1 2 10 i equals 1 3 5 7 9 for i 20 1 1 i equals 20 19 2 1 Note This is MATLAB built in language T 3 STATEMENTS while Purpose Define a conditional loop Syntax while logical expression end Description while initiates a conditional loop which terminates when logical expression equals False The while loop is closed by end to define th
81. f i M m end The finite element mesh is plotted using the following commands clf eldraw2 Ex Ey 1 2 2 Edof grid title 2D Frame Structure pause 2 D Frame Structure 3 4 2 57 0 57 0 5 1 1 5 2 2 5 Ox 0 5 Finite element mesh 9 4 3 EXAMPLES exdl Dynamic analysis A standard procedure in dynamic analysis is eigenvalue analysis This is accom plished by the following set of commands b 1 2 3 14 La Egv eigen K M b Freq sqrt La 2 pi Note that the boundary condition matrix b only lists the degrees of freedom that are zero The results of these commands are the eigenvalues stored in La and the eigenvectors stored in Egv The corresponding frequencies in Hz are calculated and stored in the column matrix Freq Freq 6 9826 43 0756 66 5772 162 7453 230 2709 295 6136 426 2271 697 7628 877 2765 955 9809 1751 3 The eigenvectors can be plotted by entering the commands below figure 1 clf grid title The first eigenmode eldraw2 Ex Ey 2 3 1 Edb extract Edof Egv 1 FreqText num2str Freq 1 pause eldisp2 Ex Ey Edb 1 2 2 text 5 1 75 FreqText The first eigenmode cae Ma oe ee i L 2 5M i i I 2r ry I J 6 983 1 5 1 i 1H 0 5 ee ee S 0 5 0 0 5 1 1 5 2 2 5 The first eigenmode 6 98 Hz EXAMPLES 9 4 4 Dynamic analysis exdl An attractive way of displaying the eigen
82. gt A 25 0e 4 ex1 0 2 ex6 4 4 ey1 2 2 ey6 0 2 E 2 1e11 ep E Al ex2 0 2 ex3 2 4 ex4 2 4 ex5 2 2 ex7 0 2 ex8 2 4 ex9 0 2 ex10 2 4 ey2 0 0 ey3 2 2 ey4 0 0 ey5 0 2 ey7 0 2 ey8 0 2 ey9 2 0 ey10 2 0 Kel bar2e exl eyl ep Ke2 bar2e ex2 ey2 ep Ke3 bar2e ex3 ey3 ep Ke4 bar2e ex4 ey4 ep Ke5 bar2e ex5 ey5 ep Ke6 bar2e ex6 ey6 ep Ke7 bar2e ex7 ey7 ep Ke8 bar2e ex8 ey8 ep Ke9 bar2e ex9 ey9 ep Kel0 bar2e ex10 ey10 ep K assem Edof 1 K Ke1 K assem Edof 2 K Ke2 K assem Edof 3 K Ke3 K assem Edof 4 K Ke4 K assem Edof 5 K Ke5 K assem Edof 6 K Ke6 K assem Edof 7 K Ke7 K assem Edof 8 K Ke8 K assem Edof 9 K Ke9 K assem Edof 10 K Ke10 The system of equations is solved considering the boundary conditions in bc gt gt bc 1 0 2 0 3 0 4 0 gt gt a solveq K f bc a oOoo0oo0ooO 0024 0045 0016 0042 0030 0107 0017 0113 9 3 13 EXAMPLES exs4 Static analysis The displacement at the point of loading is 1 7 1073 m in the x direction and 11 3 1073 m in the y direction Normal forces are evaluated from element displacements These are obtained from the global displacements a using the function extract The normal forces are evaluated using the function bar2s gt gt edi extract Edof 1 a gt gt ed3 extract Edof 3 a gt gt ed5 extract Edof 5
83. he element topology is defined by the topology matrix gt gt Edof 1 1 2 3 4 5 6 2 4 5 6 7 8 9 3 7 8 9 10 11 12 The system matrices i e the stiffness matrix K and the load vector f are defined by gt gt K zeros 12 f zeros 12 1 f 5 10000 The element property vector ep the element coordinate vectors ex and ey and the element stiffness matrix Ke are generated Note that the same coordinate vectors are applicable for all elements because they are identical 9 3 9 EXAMPLES exs3 Static analysis gt gt E 2 1e11 A 45 3e 4 T 2510e 8 ep LE A I gt gt ex 0 3 ey L0 0 gt gt Ke beam2e ex ey ep Ke 1 0e 08 3 1710 O 0 3 1710 0 0 0 0 0234 0 0351 0 0 0234 0 0351 0 0 0351 0 0703 0 0 0351 0 0351 3 1710 0 0 3 1710 0 0 0 0 0234 0 0351 0 0 0234 0 0351 0 0 0351 0 0351 0 0 0351 0 0703 Based on the topology information the global stiffness matrix can be generated by assembling the element stiffness matrices gt gt K assem Edof K Ke Finally the solution can be calculated by defining the boundary conditions in bc and solving the system of equations Displacements a and reaction forces Q are computed by the function solveq gt gt bc 1 0 2 0 11 0 a Q solveq K f bc The section forces es are calculated from element displacements Ed gt gt Ed extract Edof a gt gt esl beam2s ex ey ep Ed 1 gt gt es2 beam2s ex ey ep Ed 2 gt gt es3 beam2s ex
84. inate matrix eigen Solve a generalized eigenvalue problem extract Extract values from a global vector insert Assemble element internal force vector solveq Solve a system of equations statcon Perform static condensation 6 2 1 SYSTEM assem Static system functions Purpose Assemble element matrices l J ky kp i J k ka ky eee eee k ke ky k urunk i ki k o j Raa k k mea j kK i dof r J z dof nn K Syntax K assem edof K Ke K f assem edof K Ke f fe Description assem adds the element stiffness matrix K stored in Ke to the structure stiffness matrix K stored in K according to the topology matrix edof The element topology matrix edof is defined as edof el dof dof dofnedl M global dof where the first column contains the element number and the columns 2 to ned 1 contain the corresponding global degrees of freedom ned number of element de grees of freedom In the case where the matrix K is identical for several elements assembling of these can be carried out simultaneously Each row in Edof then represents one element i e nel is the total number of considered elements ely dof do f2 o wa 43 do fned el do do df ne Edof fa fz f i one row for each element elnel dof do fz Ge do fned If fe and f are given in the function the element load vector f is also added to the global load vector f SYSTEM 6 2 2 Static syst
85. ines for the plange element follows the numerical approximation along the element boundaries i e a linear variation A finer element mesh will bring the contour lines closer to the corner point Along the symmetry line the contour lines should physically be perpendicular to the boundary This will also be improved with a finer element mesh With the MATLAB functions colormap and fill a color plot of the concentrations can be obtained colormap jet fill Ex Ey Ed axis equal EXAMPLES 9 3 32 9 4 Dynamic analysis This section concerns linear dynamic finite element calculations The examples illustrate some basic features in dynamics such as modal analysis and time stepping procedures Dynamic analysis exd1 exd2 exd3 exd4 Modal analysis of frame Transient analysis Reduced system transient analysis Time varying boundary condition Note The examples listed above are supplied as m files on the CALFEM diskette under the directory examples The example files are named according to the table 9 4 1 EXAMPLES exdl Dynamic analysis Purpose Set up the finite element model and perform eigenvalue analysis for a simple frame structure Description Consider the two dimensional frame shown below A vertical beam is fixed at its lower end and connected to a horizontal beam at its upper end The horizontal beam is simply supported at the right end The length of the vertical
86. ing in workspace in a binary file named file name mat save filename variables writes named variables separated by blanks in a binary file named filename mat save filename variables ascii writes named variables in an ASCII file named filename Note This is a MATLAB built in function For more information about the save function type help save GENERAL PURPOSE 2 10 type Purpose List file Syntax type filename Description type filename lists the specified file Use path names in the usual way for your operating system If a filename extension is not given m is added by default This makes it convenient to list the contents of m files on the screen Note This is a MATLAB built in function For more information about the type function type help type 2 11 GENERAL PURPOSE what Purpose Directory listing of m files mat files and mex files Syntax what what dirname Description what lists the m files mat files and mex files in the current directory what dirname lists the files in directory dirname in the MATLAB search path The syntax of the path depends on your operating system Note This is a MATLAB built in function For more information about the what function type help what GENERAL PURPOSE ZAI who whos Purpose List directory of variables in matrix bank Syntax who whos Description who lists the variables currently in memory whos lists the curr
87. ions of the basic equations au dO dtu EA q 0 El q 0 EI amt dg aa The equations are valid if gj is not more than a linear function of z From these equations the displacements along the beam element are obtained as the sum of the homogeneous and the particular solutions qz 0 u Z u z ur u z where azb 2 T IBAS OL m Up T T252 fame T u NC Ga u v Z w ee qgLi a O z 24E1I L 2GA ks L p T Lz oF _ QgL Z T T eS T DBT z T and g 1z0 0 0 0 EI N 10 01 x a Gag 0 0 0 1 2x 3 x 2a 1 0 0 0 0 0 0 01 0 0 0 Uy ga ee 0 6a ee gt dl L000 0 E Or 302k be L ug 0 0 0 1 2L 3 L 2a The transformation matrix G and nodal displacements a are described in beam2e Note that the transpose of a is stored in ed Finally the section forces are obtained from dit d d0 VU N EA GAk 90 M EI dz ae er dz ELEMENT 5 6 10 Two dimensional beam element beam2w Purpose Compute element stiffness matrix for a two dimensional beam element on elastic foundation y lt lt k k x Syntax Ke beam2w ex ey ep Ke fe beam2w ex ey ep eq Description beam2w provides the global element stiffness matrix Ke for a two dimensional beam element on elastic foundation The input variables ex and ey are described in beam2e and ep E A I ka k contains the modulus of elasticity the cross section area A the moment of inertia J
88. is computed according to f Blo tdA A where the matrices B and o are defined in plani4e and plani4s respectively Evaluation of the integral is done by Gauss integration 5 9 25 ELEMENT plani8e Two dimensional solid elements Purpose Compute element matrices for an 8 node isoparametric element in plane strain or plane stress 7 3 4 8 6 y u 1 5 2 x Y Ui Syntax Ke plani8e ex ey ep D Ke fe plani8e ex ey ep D eq Description plani8e provides an element stiffness matrix Ke and an element load vector fe for an 8 node isoparametric element in plane strain or plane stress The element nodal coordinates 71 y1 2 etc are supplied to the function by ex and ey The type of analysis ptype the element thickness t and the number of Gauss points n are supplied by ep ptype 1 plane stress n x n integration points ptype 2 plane strain n 1 2 3 The material properties are supplied by the constitutive matrix D Any arbitrary D matrix with dimensions from 3 x 3 to 6 x 6 may be given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 ex 21 2 g ep ptype t n ey y Yo Ys ptyp Dy Di2 Di3 D D s Di6 Do Dog Do3 D Do5 D26 Dy D D E 11 12 13 B D31 D2 Ds 3 D34 D5 D36 D Da Do Do3 or D Day D32 D43 Daa Das D46 D3 D32 D33 Ds Ds2 Ds3 Ds4 Dss Dse Dei Dez Des Dea
89. is obtained by the function extract as ed a u uz ug If body forces are applied to the element the variable eq must be included The output variables es o7 Oxz Oyy oz Oxy ozz yz et et Exa Eyy Ezz Yay z2 ye contain the stress and strain components The size of es and et follows the size of D Note that for plane stress 0 and for plane strain o 0 Theory By assembling triangular elements as described in plange a system of equations con taining 10 degrees of freedom is obtained From this system of equations the two unknown displacements at the center of the element are computed Then according to the description in plants the strain and stress components in each of the four trian gular elements are produced Finally the quadrilateral element strains and stresses are computed as area weighted mean values from the values of the four triangular elements If uniformly distributed loads are applied on the element the element load vector eq is needed for the calculations 5 5 11 ELEMENT planre Two dimensional solid elements Purpose Compute element matrices for a rectangular Melosh element in plane strain or plane stress uz Us Xp 4 X33 A Ug uz gt Oo Xo Syntax Ke planre ex ey ep D Ke fe planre ex ey ep D eq Description planre provides an element stiffness matrix Ke and an element load vector fe for a rectangular Mel
90. ix sparsity structure Syntax spy S Description spy S plots the sparsity structure of any matrix S S is usually a sparse matrix but the function also accepts full matrices and the nonzero matrix elements are plotted Note This is a MATLAB built in function For more information about the spy function type help spy MATRIX 3 18 sqrt Purpose Square root Syntax B sqrt A Description B sqrt A computes the square root of the elements in matrix A and stores the result in matrix B Note This is a MATLAB built in function For more information about the sqrt function type help sqrt 3 19 MATRIX sum Purpose Sum of the elements of a matrix Syntax b sum A Description For a vector a the statement b sum a results in a scalar a containing the sum of all elements of a For a matrix A the statement b sum A returns a row vector b containing the sum of the elements found in each column vector of A The sum of all elements of a matrix is determined by c sum sum A Note This is a MATLAB built in function For more information about the sum function type help sum MATRIX 3 20 Zeros Purpose Generate a zero matrix Syntax A zeros m n Description A zeros m n results in an m times n matrix A of zeros Note This isa MATLAB built in function For more information about the zeros function type help zeros MATRIX 4 Material functions
91. l coordinate system located at the center of the element the element shape functions Nf N are obtained as Nf a 2 0 4 4ab a x b y Bi 1 2 2 2 2 N 55 2 v a 2 1 gg 2 where 1 1 a z3 x and b 383 y 5 0 17 ELEMENT plantce Two dimensional solid elements The matrix B is obtained as o Ox A 7 o B VN where V 0 Oy 3 2 Oy Ox Evaluation of the integrals for the Turner Clough element can be done either ana lytically or numerically by use of a 2 x 2 point Gauss integration The element load vector ff yields 8 Km 8 Ka fF abt t 8 So Ss SS ae SO ge a k ELEMENT 5 9 18 Two dimensional solid elements plantcs Purpose Compute stresses and strains in a Turner Clough element in plane strain or plane stress 4 Us 4 Ue uz Us e gt e gt Oyy pe Dan SS Eo y lbs Uy Oxy Oxx u u gt e Oxy Oyy x Syntax es et plantcs ex ey ep ed Description plantcs computes the stresses es and the strains et in a rectangular Turner Clough ele ment in plane strain or plane stress The stress and strain components are computed at the center of the element The input variables ex ey and ep are defined in plantce The vector ed contains the nodal displacements a of the element and is obtained by the function extract as ed a
92. lement defined in row 7 of Edof and nel is the total number of considered elements SYSTEM 6 2 6 Static system functions extract Example For the two dimensional beam element the extract function will extract six nodal displacements for each element given in Edof and create a matrix Ed of size nel x 6 Uy U2 UZ U4 U5 Ug Ed Uz U2 hi iG K k Uy U2 U3 U4 U5 Ug 6 2 7 SYSTEM insert Static system functions Purpose Assemble internal element forces in a global force vector fi h A latest i fri t dof z j dofi L J f Syntax f insert edof f ef Description insert adds the internal element load vector ff stored in ef to the global internal force vector f stored in f according to the topology matrix edof The function is for use in nonlinear analysis The element topology matrix edof is defined in assem The vector f is the global internal force vector and the vector ef is the internal element force vector computed from the element stresses see for example plani4f SYSTEM 6 2 8 Static system functions solveq Purpose Solve equation system Syntax a solveq K f a solveq K f bc a Q solveq K f bc Description solveq solves the equation system Ka f where K is a matrix and a and f are vectors The matrix K and the vector f must be predefined The solution of the system of equations is stored in a vector a which is created by the function If some values of
93. ls for the triangular element yields K CT BD BC tA g Staal where the element area A is determined as 1 ELEMENT 5 4 4 Two dimensional heat flow elements flw2ts Purpose Compute heat flux and temperature gradients in a triangular heat flow element Syntax es et flw2ts ex ey D ed Description flw2ts computes the heat flux vector es and the temperature gradient et or corre sponding quantities in a triangular heat flow element The input variables ex ey and the matrix D are defined in flw2te The vector ed contains the nodal temperatures af of the element and is obtained by the function extract as ed af T To T3 The output variables es q q dy OT OT et VT VE Ox Oy contain the components of the heat flux and the temperature gradient computed in the directions of the coordinate axis Theory The temperature gradient and the heat flux are computed according to Vr BC ta q DVT where the matrices D B and C are described in flw2te Note that both the tem perature gradient and the heat flux are constant in the element 5 4 5 ELEMENT flw2qe Two dimensional heat flow elements Purpose Compute element stiffness matrix for a quadrilateral heat flow element T T X3Y3 X4y4 x ay 2 T Syntax Ke flw2qe ex ey ep D Ke fe flw2qe ex ey ep D eq Description flw2qe provides the element stiffness conductivity matrix Ke
94. matrix K For most of the elements an element load vector f can also be computed These functions are identified by their last letter e Using the function assem the element stiffness matrices and element load vectors are assembled into a global stiffness matrix K and a load vector f Unknown nodal values of temperatures or displacements a are computed by solving the system of equations Ka f using the function solveq A vector of nodal values of temperatures or displacements for a specific element are formed by the function extract When the element nodal values have been computed the element flux or element stresses can be calculated using functions specific to the element type concerned These functions are identified by their last letter s For some elements a function for computing the internal force vector is also available These functions are identified by their last letter f 51 1 ELEMENT 5 2 Spring element The spring element shown below can be used for the analysis of one dimensional spring systems and for a variety of analogous physical problems Uz Quantities corresponding to the variables of the spring are listed in Table 1 Problem type Spring Nodal dis Element Spring stiffness placement force force Spring k u P N EA Bar T u P N Thermal conduction ua T Q Q ae 1 7 Electrical circuit R U I I kA Groundwater flow T o Q Q D 2 Pipe network Pr p Q Q
95. matrix C stored in the variables Ke Me and Ce respectively are computed according to K G K G M GMG C G C G where G and K are described in beam2e The matrix M is given by 140 0 0 7 0 0 0 156 222 0 54 SBE ye mE 0 2L 4P 0 13L 3p 420 70 0 0 140 0 0 0 54 182 0 156 22L 0 13L 3L 0 22L 41 and the matrix C is computed by combining K and M Co M aik 5 6 21 ELEMENT beam2ds Two dimensional beam element Purpose Compute section forces for a two dimensional beam element in dynamic analysis V2 X M lt N J Syntax es beam2ds ex ey ep ed ev ea Description beam2ds computes the section forces at the ends of the dynamic beam element beam2d The input variables ex ey and ep are defined in beam2d The element displace ments the element velocities and the element accelerations stored in ed ev and ea respectively are obtained by the function extract The output variable MUM No Vo M contains the section forces at the ends of the beam ELEMENT 5 6 22 Two dimensional beam element beam2ds Theory The section forces at the ends of the beam are obtained from the element force vector P N Vi M N Vo M computed according to P K Ga C Gae M G a The matrices K and G are described in beam2e and the matrices M and C are described in beam2d The nodal displacements e T af u U2 uz Us Us use shown in beam2d al
96. ment k u A uz gt VW Syntax Ke springle ep Description springle provides the element stiffness matrix Ke for a spring element The input variable ep K supplies the spring stiffness k or the analog quantity defined in Table 1 Theory The element stiffness matrix K stored in Ke is computed according to te SE k where k is defined by ep ELEMENT 52 4 Spring element spring1s Purpose Compute spring force in a spring element VY Syntax es springls ep ed Description springls computes the spring force es in a spring element The input variable ep is defined in springle and the element nodal displacements ed are obtained by the function extract The output variable es N contains the spring force N or the analog quantity Theory The spring force N or analog quantity is computed according to N k w u 5 2 5 ELEMENT Spring element spring1s 5 3 Bar elements Bar elements are available for one two and three dimensional analysis dimensional element see the spring element For the one Bar elements U4 U3 Uy u bar2e bar3e bar2g Two dimensional bar functions bar2e Compute element matrix bar2g Compute element matrix for geometric nonlinear element bar2s Compute normal force Three dimensional bar functions bar3e Compute element matrix bar3s Compute normal force 5 3 1 ELEMENT
97. mestep is chosen as At 0 001 seconds and the integration is performed for T 1 0 second At every 0 1 second the deformed shape of the whole structure shall be displayed The first eigenmode a ee a HS yeaa cae l 2 5 i 2 I 6 983 gt 1 5 1 1H 0 5 SO se ees ee ee ae 0 5 0 0 5 1 15 2 2 5 Time history of the impact load The load is generated using the gfunc function The time integration is performed by the step2 function Because there is no damping present the C matrix is entered as dt 0 005 T 1 hoem the load aa aa a a G 0 0 0 15 1 0 250 T 0 t g gfunc G dt f zeros 15 length g f 4 1000 g boundary condition initial condition bc 1 0 2 0 3 0 14 0 d0 zeros 15 1 v0 zeros 15 1 output parameters ntimes 0 1 0 1 1 nhist 4 11 4 time integration parameters ip dt T 0 25 0 5 10 2 ntimes nhist time integration k sparse K m sparse M EXAMPLES 9 4 6 Dynamic analysis2 exd2 Dsnap D V A step2 k m d0 v0 ip f bc The requested time history plots are generated by the following commands figure 1 plot t D 1 t D 2 grid xlabel time sec ylabel displacement m title Displacement time for the 4th and 11th gt degree of f
98. metric heat flow element The element nodal coordinates x71 Y1 21 2 etc are supplied to the function by ex ey and ez The number of Gauss points n nx nxn integration points n 1 2 3 are supplied to the function by ep and the thermal conductivities or corresponding quantities krs Kkxy etc are supplied by D ex T T2 T3 Tg kra kry kzz ey y1 Yo ya Ys ep n D kys kyy kyz ez 21 22 23 2g k kzy kzz If the scalar variable eq is given in the function the element load vector fe is com puted using eq Q where Q is the heat supply per unit volume Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K B DB av V f NTQav V ELEMENT 5 4 14 Three dimensional heat flow elements flw3i8e with the constitutive matrix D defined by D The evaluation of the integrals for the 3D isoparametric 8 node element is based on a temperature approximation T 7 expressed in a local coordinates system in terms of the nodal variables T to Tg as T n Neat where N N NE NE NE a amp Ti Th Ts Ta The element shape functions are given by n ipini N arona Ng 04 904 H 0 9 NP 0 9a m 0 9 N s0 9 MA 9 MEZUA MUHE NS TETE EE N S ELEM The B matrix is given by a a Ox de e e o e T 1 0 BS VN N J N Oy cae On a Oz aC where J is
99. modes is shown in the figure below The result is accomplished by translating the different eigenmodes in the x direction see the Ext matrix defined below and in the y direction see the Eyt matrix clf axis equal hold on axis off magnfac 0 5 title The first eight eigenmodes Hz for i1 1 4 Ext Ex i 1 3 eldraw2 Ext Ey 2 3 1 Edb extract Edof Egv i eldisp2 Ext Ey Edb 1 2 2 magnfac FreqText num2str Freq i text 3 i 1 5 1 5 FreqText end Eyt Ey 4 for i 5 8 Ext Ex i 5 3 eldraw2 Ext Eyt 2 3 1 Edb extract Edof Egv i eldisp2 Ext Eyt Edb 1 2 2 magnfac FreqText num2str Freq i text 3 i 5 5 2 5 FreqText end The first eight eigenmodes Hz 6 983 66 58 162 7 xX I I 230 3 b 295 6 426 2 697 8 I The first eight eigenmodes Frequencies are given in Hz 9 4 5 EXAMPLES exd2 Dynamic analysis Purpose The frame structure defined in exd1 is exposed in this example to a transient load The structural response is determined by a time stepping procedure Description The structure is exposed to a transient load impacting on the center of the vertical beam in horizontal direction i e at the 4th degree of freedom The time history of the load is shown below The result shall be displayed as time history plots of the 4th degree of freedom and the 11th degree of freedom At time t 0 the frame is at rest The ti
100. mples Section forces or element displacements can easily be plotted The bending moment M along the beam is plotted by gt gt plot eci es 3 ELEMENT 5 6 6 Two dimensional Timoshenko beam element beam2t Purpose Compute element stiffness matrix for a two dimensional Timoshenko beam element X Y E G A L k a 3 Xpy Syntax Ke beam2t ex ey ep Ke fe beam2t ex ey ep eq Description beam2t provides the global element stiffness matrix Ke for a two dimensional Timo shenko beam element The input variables ex z 2 EGAI k ey y y2 a supply the element nodal coordinates 1 Y1 2 and yo the modulus of elasticity E the shear modulus G the cross section area A the moment of inertia J and the shear correction factor ks The element load vector fe can also be computed if uniformly distributed loads are applied to the element The optional input variable eq a as then contains the distributed loads per unit length qz and qz 5 6 7 ELEMENT beam2t Two dimensional Timoshenko beam element Theory The element stiffness matrix K stored in Ke is computed according to K G K G where G is described in beam2e and K is given by EA EA i 0 0 0 0 0 12EI 6EI 0 _ BEI 6EI L 1 u LA p L LL 1 p 0 6EI 4EI 1 5 0 SEI 2EI 1 5 Ke L Ltn L L p Z4 0 0 ft 0 0 0 IEI __ 651 0
101. n shape axis equal changes the current axis box size so that equal tick mark increments on the x and y axes are equal in size This makes plot sin x cos x look like a circle instead of an oval axis normal restores the current axis box to full size and removes any restrictions on the scaling of the units This undoes the effects of axis square and axis equal axis off turns off all axis labeling and tick marks axis on turns axis labeling and tick marks back on Note This is a MATLAB built in function For more information about the axis function type help axis GRAPHICS 8 2 clf Purpose Clear current figure graph window Syntax clf Description clf deletes all objects axes from the current figure Note This is a MATLAB built in function For more information about the clf function type help clf 8 3 GRAPHICS eldraw2 Purpose Draw the undeformed mesh for a two dimensional structure Syntax eldraw2 Ex Ey eldraw2 Ex Ey plotpar eldraw2 Ex Ey plotpar elnum Description eldraw2 displays the undeformed mesh for a two dimensional structure Input variables are the coordinate matrices Ex and Ey formed by the function co ordxtr The variable plotpar sets plot parameters for linetype linecolor and node marker plotpar linetype linecolor nodemark linetype 1 solid line linecolor 1 black 2 dashed line 2 blue 3 dotted line 3 magenta 4 red nodemark 1 circle 2
102. nal that can be identified by the Fourier transform corresponds to half the sampling frequency The sampling frequency is equal to 1 dt where dt is the time increment of the time signal The complex Fourier coefficients p k are stored in the vector p and are computed according to j 1 k 1 plk 2 e J wn Me 1 j where wy eT 2nt N Note This is a MATLAB built in function SYSTEM 6 3 4 Dynamic system functions freqresp Purpose Compute frequency response of a known discrete time response Syntax Freq Resp freqresp D dt Description freqresp computes the frequency response of a discrete dynamic system D is the time history function and dt is the sampling time increment i e the time increment used in the time integration procedure Resp contains the computed response as a function of frequency Freq contains the corresponding frequencies Example The result can be visualized by gt gt plot Freq Resp gt gt xlabel frequency Hz or gt gt semilogy Freq Resp gt gt xlabel frequency Hz The dimension of Resp is the same as that of the original time history function Note The time history function of a discrete system computed by direct integration behaves often in an unstructured manner The reason for this is that the time history is a mixture of several participating eigenmodes at different eigenfrequencies By using a Fourier transform how
103. ne strain or plane stress This element can only be used if the material is isotropic and if the element edges are parallel to the coordinate axis The element nodal coordinates x1 y1 and 3 y3 are supplied to the function by ex and ey The state of stress ptype the element thickness t and the material properties E and v are supplied by ep For plane stress ptype 1 and for plane strain ptype 2 ex z zz ep pype t E v ey y ys p ptyp If uniformly distributed loads are applied to the element the element load vector fe is computed The input variable be a containing loads per unit volume by and b is then given ELEMENT 5 9 16 Two dimensional solid elements plantce Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT DB dA f NT bia where the constitutive matrix D is described in hooke see Section 4 and the body force vector b is defined by eq The evaluation of the integrals for the Turner Clough element is based on a dis placement field u x y built up of a bilinear displacement approximation superposed by bubble functions in order to create a linear stress field over the element The displacement field is expressed in terms of the nodal variables u1 us Ug as u z y N a where ui ug Ne Ne NE INE NS NS NE N NS a u2 saa FARNE ONE NG ING Ne NS VE Ne gE Ug With a loca
104. near system of equations Ka f where K is the global stiffness matrix and f is the global load vector Often used static system functions are assem and solveq The function assem assembles the global stiffness matrix and solveq computes the global displacement vector a considering the boundary conditions It should be noted that K f and a also represent analogous quantities in systems others than structural mechanical systems For example in a heat flow problem K represents the conductivity matrix f the heat flow and a the temperature Dynamic system functions are related to different aspects of linear dynamic systems of coupled ordinary differential equations according to CdiKd f for first order systems and D Md Cd Kd f for second order systems First order systems occur typically in transient heat conduction and second order systems occur in structural dynamics 6 1 1 SYSTEM 6 2 Static system functions The group of static system functions comprises functions for setting up and solving the global system of equations It also contains a function for eigenvalue analysis a function for static condensation a function for extraction of element displacements from the global displacement vector and a function for extraction of element coordinates The following functions are available for static analysis Static system functions assem Assemble element matrices coordxtr Extract element coordinates from a global coord
105. ntf Two dimensional solid elements Purpose Compute internal element force vector in a triangular element in plane strain or plane stress Syntax ef plantf ex ey ep es Description plantf computes the internal element forces ef in a triangular element in plane strain or plane stress The input variables ex ey and ep are defined in plante and the input variable es is defined in plants The output variable ef f fa fiz Shs fis contains the components of the internal force vector Theory The internal force vector is computed according to f J Bla t dA A where the matrices B and are defined in plante and plants respectively Evaluation of the integral for the triangular element yields f BC otA ELEMENT 5 0 8 Two dimensional solid elements plange Purpose Compute element matrices for a quadrilateral element in plane strain or plane stress U6 Syntax Ke planqe ex ey ep D Ke fe planqe ex ey ep D eq Description plange provides an element stiffness matrix Ke and an element load vector fe for a quadrilateral element in plane strain or plane stress The element nodal coordinates 1 Y1 2 etc are supplied to the function by ex and ey The type of analysis ptype and the element thickness t are supplied by ep ptype 1 plane stress ptype 2 plane strain and the material properties are supplied by the constitutive matrix D Any arbitrary D matrix with dimensions fr
106. nts flw2te Purpose Compute element stiffness matrix for a triangular heat flow element T x py 3 y T T Xo xpy Syntax Ke flw2te ex ey ep D Ke fe flw2te ex ey ep D eq Description flw2te provides the element stiffness conductivity matrix Ke and the element load vector fe for a triangular heat flow element The element nodal coordinates 71 Y1 2 etc are supplied to the function by ex and ey the element thickness t is supplied by ep and the thermal conductivities or corresponding quantities kss kz etc are supplied by D ex z1 T2 T3 kzz ky t D eyv 41 Yo yz P kys kyy If the scalar variable eq is given in the function the element load vector fe is com puted using eq Q where Q is the heat supply per unit volume Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to Ke C71 fB DBtdA C C717 N Qtaa with the constitutive matrix D defined by D The evaluation of the integrals for the triangular element is based on the linear temperature approximation T x y and is expressed in terms of the nodal variables Th To and Ts as T x y N a N C a 5 4 3 ELEMENT flw2te Two dimensional heat flow elements where 7 la y T Ne zxz y C 1 r2 2 ao T 1 23 Ys T3 and hence it follows that oO 010 Ox BN panel Oy Evaluation of the integra
107. o dimensional diffusion Description i 3 3 c 10 kg m 0 1m E E it D gt 7 8 10 ju 12 5 6 c 0 7 8 9 3 4 2 gt L6 1 2 ie 3 c 0 O lm x Description Consider a filter paper of square shape Three sides are in contact with pure water and the fourth side is in contact with a solution of concentration c kg m The length of each side is 0 100 m Using symmetry only half of the paper has to be analyzed The paper and the corresponding finite element mesh are shown The following boundary conditions are applied c 0 y clx 0 c 0 1 y 0 c x 0 1 107 The element topology is defined by the topology matrix gt gt Edof 1 1 2 5 4 2 2 3 6 5 3 4 5 8 7 4 5 6 9 8 5 7 8 11 10 6 8 912 11 7 10 11 14 13 8 11 12 15 14 9 3 27 EXAMPLES exs7 Static analysis The system matrices i e the stiffness matrix K and the load vector f are defined by gt gt K zeros 15 f zeros 15 1 Because of the same geometry orientation and constitutive matrix for all elements only one element stiffness matrix Ke has to be computed This is done by the function flw2qe gt gt ep 1 D 1 0 0 1 gt gt ex 0 0 025 0 025 0 ey 0 0 0 025 0 025 gt gt Ke flw2qe ex ey ep D gt gt Ke 0 7500 0 2500 0 2500 0 2500 0 2500 0 7500 0 2500 0 2500 0 2500 0 2500 0 7500 0 2500 0 2500 0 2500 0 2500 0 7500 Based on the topology information the global stiffness m
108. of Gauss points n n xn integration points n 1 2 3 are supplied to the function by ep and the thermal conductivities or corresponding quantities kss ky etc are supplied by D ex t1 T2 T3 Tg kzz aa ep t n D y ey y1 Yo ys Ys P Le kyy If the scalar variable eq is given in the function the vector fe is computed using eq Q where Q is the heat supply per unit volume 5 4 11 ELEMENT flw2i8e Two dimensional heat flow elements Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT DB dA f NT Quad with the constitutive matrix D defined by D The evaluation of the integrals for the 2D isoparametric 8 node element is based on a temperature approximation T 7 expressed in a local coordinates system in terms of the nodal variables T to Tg as T 7 Nfae where N N NE NE NE a Ti Th Tz Ta The element shape functions are given by NE IQA N E n 1 1 Ny Gll 1 n 1 n Ne 5 1 6 0 1 1 N3 70 En E n N EELA 1 1 Ni 30 6 m n Np The Bmatrix is given by o o e_ e Ox e T 1 o e Be VNS a N J a N Oy On where J is the Jacobian matrix r de ya 8 On Oy Oy Of On Evaluation of the integrals is done by Gauss integration ELEMENT 5 4 12 Two dimensional heat flow elements flw2i8s Purpose
109. of elasticity EF the cross section area A and the moment of inertia J The element load vector fe can also be computed if uniformly distributed loads are applied to the element The optional input variable eq a qs then contains the distributed loads per unit length qz and qj ELEMENT 5 6 2 Two dimensional beam element beam2e 1 x Theory The element stiffness matrix K stored in Ke is computed according to K G K G where EA Q 0 0 0 EI EI EI EI Re 0 2 e 4 o 0 FA 9 0 L L 0 mee _ 0 ea ae 6EI 2EI 6EI 4EI OR zr 8 p a Neg Ny 0O 0 0 0 0 0 1 O0 0 o gt 0 0 0 n rye 0 0 0 0 Nag Nyy 0 0 0 0 0 0 1 The transformation matrix G contains the direction cosines T2 T maen _ Y2 7y z liry L d L Nee Nyg where the length L V x2 21 y2 m1 The element load vector ff stored in fe is computed according to fe GTF 5 6 3 ELEMENT beam2e Two dimensional beam element where ELEMENT 5 6 4 Two dimensional beam element beam2s Purpose Compute section forces in a two dimensional beam element Va M N A V sik a M NN 7 pe C N V Syntax es beam2s ex ey ep ed es beam2s ex ey ep ed eq es edi eci beam2s ex ey ep ed eq n Description beam2s computes the section forces and displacements in local directions along the beam element beam2e The input variables ex ey ep and eq are defin
110. ollowing functions are described in this chapter Vector and matrix operations Special characters Special characters Create vectors and do matrix subscripting x Matrix arithmetic abs Absolute value det Matrix determinant diag Diagonal matrices and diagonals of a matrix inv Matrix inverse length Vector length max Maximum element s of a matrix min Minimum element s of a matrix ones Generate a matrix of all ones red Reduce the size of a square matrix size Matrix dimensions sqrt Square root sum Sum of the elements of a matrix zeros Generate a zero matrix Sparse matrix handling full Convert sparse matrix to full matrix sparse Create sparse matrix spy Visualize sparsity structure MATRIX lO 53 Purpose Special characters Syntax LOH a5 Description Brackets are used to form vectors and matrices Parentheses are used to indicate precedence in arithmetic expressions and to specify an element of a matrix Used in assignment statements Matrix transpose X is the transpose of X If X is complex the apostrophe sign performs complex conjugate as well Do X if only the transpose of the complex matrix is desired Decimal point 314 100 3 14 and 0 314e1 are all the same Comma Used to separate matrix subscripts and function arguments Semicolon Used inside brackets to end rows Used after an expression to suppress printing or to separate statements
111. om 3 x 3 to 6 x 6 may be given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 ex 21 T2 3 X4 ep ptype t ey y y2 Ys ya p pipe t D D D 7 u Dig Dis Ds D32 D33 Da D35 D36 D3 D32 D33 5 9 9 ELEMENT plange Two dimensional solid elements If uniformly distributed loads are applied on the element the element load vector fe is computed The input variable by el containing loads per unit volume by and b is then given Theory In computing the element matrices two more degrees of freedom are introduced The location of these two degrees of freedom is defined by the mean value of the coordinates at the corner points Four sets of element matrices are calculated using plante These matrices are then assembled and the two extra degrees of freedom are eliminated by static condensation ELEMENT 5 9 10 Two dimensional solid elements plangs Purpose Compute stresses and strains in a quadrilateral element in plane strain or plane stress Syntax es et planqs ex ey ep D ed es et planqs ex ey ep D ed eq Description plangs computes the stresses es and the strains et in a quadrilateral element in plane strain or plane stress The input variables ex ey ep D and eq are defined in plange The vector ed contains the nodal displacements a of the element and
112. on 5 0 37 ELEMENT 5 6 Beam elements Beam elements are available for two and three dimensional linear static analysis Two dimensional beam elements for nonlinear geometric and dynamic analysis are also available Beam elements beam2e beam2t beam2w beam3e beam2g beam2d 2D beam elements beam2e Compute element matrices beam2s Compute section forces beam2t Compute element matrices for Timoshenko beam element beam2ts Compute section forces for Timoshenko beam element beam2w Compute element matrices for beam element on elastic foundation beam2ws Compute section forces for beam element on elastic foundation beam2g Compute element matrices with respect to geometric nonlinearity beam2gs Compute section forces for geometric nonlinear beam element beam2d Compute element matrices dynamic analysis beam2ds Compute section forces dynamic analysis 3D beam elements beam3e Compute element matrices beam3s Compute section forces 5 6 1 ELEMENT beam2e Two dimensional beam element Purpose Compute element stiffness matrix for a two dimensional beam element Xpy Syntax Ke beam2e ex ey ep Ke fe beam2e ex ey ep eq Description beam2e provides the global element stiffness matrix Ke for a two dimensional beam element The input variables ex z x ey y y pi aed supply the element nodal coordinates 71 y1 2 and y2 the modulus
113. ons platre Compute element matrices platrs Compute section forces 5 7 1 ELEMENT platre Plate element Purpose Compute element stiffness matrix for a rectangular plate element Uj Ug u u Xa meti 2 i Uy X3 3 y u3 Us Uy i CY P lt u j E2 ea tA a Syntax Ke platre ex ey ep D Ke fe platre ex ey ep D eq Description platre provides an element stiffness matrix Ke and an element load vector fe for a rectangular plate element This element can only be used if the element edges are parallel to the coordinate axis The element nodal coordinates x1 y1 2 etc are supplied to the function by ex and ey the element thickness t by ep and the material properties by the constitutive matrix D Any arbitrary D matrix with dimensions 3 x 3 and valid for plane stress may be given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 7 Dy Diy Dig ex z T2 T3 z4 ep t D Da D Dog ey Yo y3 ya Ds D32 D33 If a uniformly distributed load is applied to the element the element load vector fe is computed The input variable eq q then contains the load q per unit area in the z direction ELEMENT 5 7 2 Plate element platre Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to Ke C71
114. order differential equa tions of the form Md Cd Kd f x t d 0 do d 0 vo In structural mechanics problems K C and M represent the n x n stiffness damping and mass matrices respectively The initial conditions are given by the vectors d0 and vO containing initial dis placements and initial velocities The time integration procedure is governed by the parameters given in the vector ip defined as ip dt Tad nsnap nhist time dof e a _ eee list of list of nsnap nhist moments dofs where dt specifies the time increment in the time stepping scheme T the total time and a and time integration constants for the Newmark family of methods see 1 The parameter nsnap denotes the number of snapshots stored in Dsnap The selected elapsed times are specified in time whereas nhist is the number of time histories stored in D V and A The selected degrees of freedom are specified in dof The following table lists frequently used values of a and a i Average acceleration trapezoidal rule a i 5 Linear acceleration a 0 6 4 Central difference The matrix f contains the time dependent load vectors If no external loads are active the matrix corresponding to f should be replaced by The matrix pdisp contains the time dependent prescribed displacement If no displacements are prescribed the matrix corresponding to pdisp should be replaced by The matrix f is o
115. osh element in plane strain or plane stress This element can only be used if the element edges are parallel to the coordinate axis The element nodal coordinates x1 y and x3 y3 are supplied to the function by ex and ey The type of analysis ptype and the element thickness t are supplied by ep ptype 1 plane stress ptype 2 plane strain and the material properties are supplied by the constitutive matrix D Any arbitrary D matrix with dimensions from 3 x 3 to 6 x 6 may be given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 ex x z3 ep ptype t ey y y E Dii Dig Di3 Di4 Dis Di6 Do Do Do3 Do Dos D25 E Da Das Da STD De Dig Da Del Dg 31 D32 D33 D51 Ds2 Ds3 Ds4 Dss Dse De Dez Des Des Des Dee ELEMENT 5 9 12 Two dimensional solid elements planre If uniformly distributed loads are applied on the element the element load vector fe is computed The input variable by containing loads per unit volume by and b is then given Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT DB dA f NT bia with the constitutive matrix D defined by D and the body force vector b defined by eq The evaluation of the integrals for the rectangular element is based on a bilinear displacement approximation u x y and is ex
116. ows and columns in Y in order to avoid the complex conjugate transpose of Y see Section 3 The time response is represented by the real part of the output from the ifft command If the transpose is used and the result is stored in a matrix X each row will represent the time response for each equation as the output from the command dyna2 See also gfunc fft ifft 6 3 3 SYSTEM fft Dynamic system functions Purpose Transform functions in time domain to frequency domain Syntax p fft g p fft g N Description fft transforms a time dependent function to the frequency domain The function to be transformed is stored in the vector g Each row in g contains the value of the function at equal time intervals The function represents a span oo lt t lt however only the values within a typical period are specified by g The fft command can be used with one or two input arguments If N is not specified the numbers of frequencies used in the transformation is equal to the the numbers of points in the time domain i e the length of the variable g and the output will be a vector of the same size containing complex values representing the frequency content of the input signal The scalar variable N can be used to specify the numbers of frequencies used in the Fourier transform The size of the output vector in this case will be equal to N It should be remembered that the highest harmonic component in the time sig
117. pressed in terms of the nodal variables U1 U2 Ug as u z y N a where ui s pE Ea ELS el Us With a local coordinate system located at the center of the element the element shape functions Nf Nf are obtained as Z 1 Ni Jab z2 y ya 7 1 N Fp lt T BY 9s r 1 N3 Jag 2 T9 2 1 Ng a ET tally 1 where 1 1 5 3 x and b z3 y 5 9 13 ELEMENT planre Two dimensional solid elements The matrix B is obtained as o Ox is S B VN where V 0 Oy 3 2 Oy Ox If a larger D matrix than 3 x 3 is used for plane stress ptype 1 the D matrix is reduced to a 3 x 3 matrix by static condensation using oz Osz Oyz 0 These stress components are connected with the rows 3 5 and 6 in the D matrix respectively If a larger D matrix than 3 x 3 is used for plane strain ptype 2 the D matrix is reduced to a 3 x 3 matrix using Vez Fz 0 This implies that a 3 x 3 D matrix is created by the rows and the columns 1 2 and 4 from the original D matrix Evaluation of the integrals for the rectangular element can be done either analytically or numerically by use of a 2 x 2 point Gauss integration The element load vector ff yields 8 lt 8 8 Km 8 Sa a a w ELEMENT 5 0 14 Two dimensional solid elements planrs Purpose Compute stresses and strains in a rectangular Melosh element in plane strain or
118. ra computes the seismic response spectrum for a known acceleration history function The computation is based on the vector a that contains an acceleration time history function defined at equal time steps The time step is specified by the variable dt The value of the damping ratio is given by the variable xi Output from the computation stored in the vector s is achieved at frequencies specified by the column vector f Example The following procedure can be used to produce a seismic response spectrum for a damping ratio 0 05 defined at 34 logarithmicly spaced frequency points The acceleration time history a has been sampled at a frequency of 50 Hz corresponding to a time increment dt 0 02 between collected points gt gt freq logspace 0 1l0g10 2 33 6 34 gt gt xi 0 05 gt gt dt 0 02 gt gt s spectra a xi dt freq The resulting spectrum can be plotted by the command gt gt loglog freq s 6 3 9 SYSTEM step1 Dynamic system functions Purpose Compute the dynamic solution to a set of first order differential equations Syntax Tsnap step1 K C d0 ip f pbound Tsnap D V step1 K C d0 ip f pbound Description stepl computes at equal time steps the solution to a set of first order differential equations of the form Cd Kd f z t d 0 do The command solves transient field problems In the case of heat conduction K and C represent the n x n conductivity
119. rains for an elasto plastic isotropic hardening von Mises material Syntax es deps st mises ptype mp est st Description mises computes updated stresses es plastic strain increments deps and states vari ables st for an elasto plastic isotropic hardening von Mises material The input variable ptype is used to define the type of analysis cf hooke The vector mp contains the material constants mp Evh where F is the modulus of elasticity v is the Poisson s ratio and h is the plastic modulus The input matrix est contains trial stresses obtained by using the elas tic material matrix D in plants or some similar s function and the input vector st contains the state parameters st yi oy eey at the beginning of the step The scalar yi states whether the material behaviour is elasto plastic yi 1 or elastic yi 0 The current yield stress is denoted by a and the effectiv plastic strain by 2 gt The output variables es and st contain updated values of es and st obtained by integration of the constitutive equations over the actual displacement step The increments of the plastic strains are stored in the vector deps If es and st contain more than one row then every row will be treated by the com mand Note It is not necessary to check whether the material behaviour is elastic or elasto plastic this test is done by the function The computation is based on an Euler Backward method i e the radial return me
120. re written in a file named m file and evaluated by writing the file name in the command window The batch oriented mode is a more flexible way of performing finite element analysis because the m file can be written in an ordinary editor This way of using CALFEM is recommended because it gives a structured organization of the functions Changes and reruns are also easily executed in the batch oriented mode A command line consists typically of functions for vector and matrix operations calls to functions in the CALFEM finite element library or commands for workspace operations An example of a command line for a matrix operation is C A B where two matrices A and B are added together and the result is stored in matrix C The matrix B is the transpose of B An example of a call to the element library is Ke barle k where the two by two element stiffness matrix K is computed for a spring element with spring stiffness k and is stored in the variable Ke The input argument is given within parentheses after the name of the function Some functions have multiple input argu ments and or multiple output arguments For example lambda X eigen K M computes the eigenvalues and eigenvectors to a pair of matrices K and M The output variables the eigenvalues stored in the vector lambda and the corresponding eigenvectors stored in the matrix X are surrounded by brackets and separated by commas The input arguments are gi
121. reedom text 0 3 0 009 solid line impact point x direction text 0 3 0 007 dashed line center horizontal beam y direction Displacement time at the 4th and 11th degree of freedom solid line impact point x direction 0 015 dashed line center horizontal beam y direction 0 01 0 005 displacement m Oo 0 005 0 01 i i i 0 0 2 0 4 0 8 1 0 6 time sec Time history at DOF 4 and DOF 11 The deformed shapes at time increment 0 1 sec are stored in Dsnap They are visu alized by the following commands figure 2 clf axis equal hold on axis off magnfac 25 title Snapshots sec magnification 25 for i 1 5 Ext Ex i 1 3 eldraw2 Ext Ey 2 3 0 Edb extract Edof Dsnap i eldisp2 Ext Ey Edb 1 2 2 magnfac Time num2str ntimes i text 3 i 1 5 1 5 Time end Eyt Ey 4 for i 6 10 Ext Ex i 6 3 eldraw2 Ext Eyt 2 3 0 9 4 7 EXAMPLES exd2 Dynamic analysis Edb extract Edof Dsnap i eldisp2 Ext Eyt Edb 1 2 2 magnfac Time num2str ntimes i text 3 i 6 5 2 5 Time end Snapshots sec magnification 25 1 0 2 0 3 0 4 0 5 X X 5 aK k Ke l 0 6 l 0 7 0 8 0 9 1 x x xX Snapshots of the deformed geometry for every 0 1 sec EXAMPLES 9 4 8 Dynamic analysis exd3 Purpose This example concern
122. rganized in the following manner SYSTEM 6 3 12 Dynamic system functions step2 time history of the load at dof time history of the load at dof time history of the load at dofn The dimension of f is number of degrees of freedom x number of timesteps 1 The matrix pdisp is organized in the following manner dof time history of the displacement dofz time history of the displacement pdisp dofm time history of the displacement The dimension of pdisp is number of dofs with prescribed displacement x number of timesteps 2 The time history functions can be generated using the command gfunc If all the values of the time histories of f or pdisp are kept constant these values need to be stated only once In this case the number of columns in f is one and in pdisp two It is highly recommended to define f as sparse a MATLAB built in function In most cases only a few degrees of freedom are affected by the exterior load and hence the matrix contains only few non zero entries The computed displacement snapshots are stored in Dsnap one column for each requested snapshot according to ip i e the dimension of Dsnap is number of dofs x nsnap The computed time histories of d d and d are stored in D V and A respectively one line for each requested degree of freedom according to ip The dimension of D is nhist x number of timesteps 1 6 3 13 SYSTEM sweep Dynamic system functions Purpose Compu
123. rom these equations the displacements along the beam element are obtained as the sum of the homogeneous and the particular solutions u z v Z E p z where azla 2EA L ui Z Ge zo _ U T L u N C G a up 5 2 see z W Z caa a Pp 24EL l Go z 2GK L and 1z000 0 00 0 0 0 0 N 001zz773000000 1000000 1zz7z00 00000 00 00 01 2 ELEMENT 5 6 28 Three dimensional beam element beam3s ia Oe 0 BOE U0 0 A 0 20 30k 02 0 0 0 0 0300 000000100000 000000000010 000000010000 uy OO Ae 0 O00 0 0 ONY all TO O 0 070 1 FG EET DBPL O 00 07 0 0 0 iiss OF 0 705 O DW E 00 0 0 00 O60 0 7 E 0000 0 0 0 1 2b3L70 0 03 0 01 IR 0 0 0 0 00 The transformation matrix G and nodal displacements af are described in beam3e Note that the transpose of a is stored in ed Finally the section forces are obtained from dit do Pw N FA W EI Vz EI di 2 dz dz3 dg dw dv T GK M El Mz El di Z dz dz Examples Section forces or element displacements can easily be plotted The bending moment M along the beam is plotted by gt gt plot eci es 5 5 6 29 ELEMENT 5 7 Plate element Only one plate element is currently available a rectangular 12 dof element The element presumes a linear elastic material which can be isotropic or anisotropic Plate elements Uy Ug u1 Ug uij uz Plate functi
124. s Purpose Compute element matrices for an 8 node isoparametric solid element Syntax Ke soli8e ex ey ez ep D Ke fe soli8e ex ey ez ep D eq Description soli8e provides an element stiffness matrix Ke and an element load vector fe for an 8 node isoparametric solid element The element nodal coordinates 71 y1 Z1 2 etc are supplied to the function by ex ey and ez and the number of Gauss points n are supplied by ep n x n integration points n 1 2 3 The material properties are supplied by the constitutive matrix D Any arbitrary D matrix with dimensions 6 x 6 may be given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 Diy Dig Dis ex E T2 g Dor ses s D eyv y1 Yo Yg ep n D l l i ez 2z1 2 Z E E e De Deo Des If different D matrices are used in the Gauss points these D matrices are stored in a global vector D For numbering of the Gauss points see eci in soli8s D pel Dis If uniformly distributed loads are applied to the element the element load vector fe is computed The input variable ELEMENT 5 9 32 Three dimensional solid elements soli8e containing loads per unit volume by by and bz is then given Theory The element stiffness matrix K and the element load vector ff stored in Ke and fe respectively are computed according to K BT D B dV V ff i NT b
125. s ex ey ep D ed Description flw2i4s computes the heat flux vector es and the temperature gradient et or corre sponding quantities in a 4 node isoparametric heat flow element The input variables ex ey ep and the matrix D are defined in flw2i4e The vector ed contains the nodal temperatures a of the element and is obtained by extract as ed a Ti To Ts T The output variables g q pons q q n on ge 308 T T Ox Oy T T Ti y x et VT Ox Oy ei si a ne x Oy contain the heat flux the temperature gradient and the coordinates of the integra tion points The index n denotes the number of integration points used within the element cf flw2i4e Theory The temperature gradient and the heat flux are computed according to VT Be a q DVT where the matrices D B and a are described in flw2i4e and where the integration points are chosen as evaluation points ELEMENT 5 4 10 Two dimensional heat flow elements flw2i8e Purpose Compute element stiffness matrix for an 8 node isoparametric heat flow element Syntax Ke flw2i8e ex ey ep D Ke fe flw2i8e ex ey ep D eq Description flw2i8e provides the element stiffness conductivity matrix Ke and the element load vector fe for an 8 node isoparametric heat flow element The element nodal coordinates x1 Y1 2 etc are supplied to the function by ex and ey The element thickness t and the number
126. s in plane stress panels and plane strain and for general three dimensional analysis In the two dimensional case there are a triangular three node element a quadrilateral four node element two rectangular four node elements and quadrilateral isoparametric four and eight node elements For three dimensional analysis there is an eight node isoparametric element The elements are able to deal with both isotropic and anisotropic materials The triangular element and the three isoparametric elements can also be used together with a nonlinear material model The material properties are specified by supplying the constitutive matrix D as an input variable to the element functions This matrix can be formed by the functions described in Section 4 5 5 1 ELEMENT Solid elements Ug Us U4 U3 uz uy plante plange 4 Ug A Us Uz Us gt gt A uy A Uy uy U3 gt 5 gt planre lantce P plani4e 7 3 4 6 8 uz 1 5 2 ty plani8e soli8e ELEMENT 5 5 2 2D solid functions plante Compute element matrices for a triangular element plants Compute stresses and strains plantf Compute internal element forces plange Compute element matrices for a quadrilateral element planqs Compute stresses and strains planre Compute element matrices for a rectangular Melosh element planrs Compute stresses and strains plantce Compute element matrices for a rec
127. s reduced system analysis for the frame structure defined in exdl Transient analysis on modal coordinates is performed for the reduced system Description In the previous example the transient analysis was based on the original finite element model Transient analysis can also be employed on some type of reduced system commonly a subset of the eigenvectors The commands below pick out the first two eigenvectors for a subsequent time integration see constant nev The result in the figure below shall be compared to the result in exd2 Displacement time at the 4th and 11th degree of freedom T T 1 1 T 7 7 solid line impact point x direction 0 015 dashed line center horizontal beam y direction 0 01 0 005 displacement m 0 005 0 01 L L 1 1 L L L L ii 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 time sec Time history at DOF 4 and DOF 11 using two eigenvectors dt 0 002 T 1 nev 2 Ba he load anes asa see ee eee G 6 0 0 15 1 0 25 0 T O t g gfunc G dt f zeros 15 length g f 4 9000 g fr sparse 1 1 nev Egv 1 nev f reduced system matrices kr sparse diag diag Egv 1 nev K xEgv 1 nev mr sparse diag diag Egv 1 nev xM xEgv 1 nev 4 initial condition drO zeros nev 1 vr0 zeros nev 1 output parameters ntimes
128. segments in the loading time function See also gfunc SYSTEM 6 3 2 Dynamic system functions dyna2f Purpose Compute the dynamic solution to a set of uncoupled second order differential equa tions Syntax Y dyna2f w2 xi f p dt Description dyna2f computes the solution to the set i Uwik wea fig t ewer of differential equations in the frequency domain The vectors w2 xi and f are the squared circular frequencies w the damping ratios and the applied forces f respectively The force vector p contains the Fourier coefficients p k formed by the command fft The solution in the frequency domain is computed at equal time increments defined by dt The result is stored in the m by n matrix Y where m is the number of equations and n is the number of frequencies resulting from the fft command The dynamic solution in the time domain is achieved by the use of the command ifft Example The dynamic solution to a set of uncoupled second order differential equations can be computed by the following sequence of commands gt gt g gfunc G dt gt gt p fft g gt gt Y dyna2f w2 xi f p dt gt gt X real ifft Y where it is assumed that the input variables G dt w2 xi and f are properly defined Note that the ifft command operates on column vectors if Y is a matrix therefore use the transpose of Y The output from the ifft command is complex Therefore use Y to transpose r
129. so define the directions of the nodal velocities he t S ho a aP a U2 UZ U4 U5 g and the nodal accelerations Lak Gach Wide Sh hls aR a z z tla tis Ue Note that the transposes of a f and a are stored in ed ev and ea respectively 5 6 23 ELEMENT beam3e Three dimensional beam element Purpose Compute element stiffness matrix for a three dimensional beam element Uj Syntax Ke beam3e ex ey ez e0 ep Ke fe beam3e ex ey ez e0 ep eq Description beam3e provides the global element stiffness matrix Ke for a three dimensional beam element The input variables ex z x ey y y eo zz Yz zz ez 21 za supply the element nodal coordinates x1 y1 etc as well as the direction of the local beam coordinate system Z y Z By giving a global vector xz yz zz parallel with the positive local Z axis of the beam the local beam coordinate system is defined The variable supplies the modulus of elasticity the shear modulus G the cross section area A the moment of inertia with respect to the y axis I the moment of inertia with respect to the Z axis and St Venant torsional stiffness K ELEMENT 5 6 24 Three dimensional beam element beam3e The element load vector fe can also be computed if uniformly distributed loads are applied to the element The optional input variable eq qz Wy qz qa then contains the distributed loads
130. star 0 no mark Default is solid black lines with circles at nodes Element numbers can be displayed at the center of the element if a column vector elnum with the element numbers is supplied This column vector can be derived from the element topology matrix Edof elnum Edof 1 i e the first column of the topology matrix Limitations Supported elements are bar beam triangular three node and quadrilateral four node elements GRAPHICS 8 4 eldisp2 Purpose Draw the deformed mesh for a two dimensional structure Syntax magnfac eldisp2 Ex Ey Ed magnfac eldisp2 Ex Ey Ed plotpar eldisp2 Ex Ey Ed plotpar magnfac Description eldisp2 displays the deformed mesh for a two dimensional structure Input variables are the coordinate matrices Ex and Ey formed by the function co ordxtr and the element displacements Ed formed by the function extract The variable plotpar sets plot parameters for linetype linecolor and node marker plotpar linetype linecolor nodemark linetype 1 solid line linecolor 1 black 2 dashed line 2 blue 3 dotted line 3 magenta 4 red nodemark 1 circle 2 star 0 no mark Default is dashed black lines with circles at nodes The magnification factor magnfac is a scalar that magnifies the element displacements for visibility The magnification factor is set automatically if it is omitted in the input list Limitations Supported elements are bar beam triangular thr
131. strain ptype 2 the D matrix is reduced to a 3 x 3 matrix using Ve Aiz 0 This implies that a 3 x 3 D matrix is created by the rows and the columns 1 2 and 4 from the original D matrix Evaluation of the integrals for the triangular element yields K C B DBC tA A f ib bj bs Bye tbl where the element area A is determined as 1 1 T Y A det 1 2 2 2 1 z3 Ys ELEMENT 5 9 6 Two dimensional solid elements plants Purpose Compute stresses and strains in a triangular element in plane strain or plane stress Ug Syntax es et plants ex ey ep D ed Description plants computes the stresses es and the strains et in a triangular element in plane strain or plane stress The input variables ex ey ep and D are defined in plante The vector ed contains the nodal displacements a of the element and is obtained by the function extract as ed a u UZ ug The output variables es o7 Oza Oyy oz Oxy ozz yz r et E Exe Eyy Ezz Yay Yee Duz contain the stress and strain components The size of es and et follows the size of D Note that for plane stress 0 and for plane strain o 0 Theory The strains and stresses are computed according to e BC a o De where the matrices D B C and a are described in plante Note that both the strains and the stresses are constant in the element 5 5 7 ELEMENT pla
132. t based graphics Mesh plots displacements section forces flows iso lines and principal stresses can be displayed The functions are divided into two dimensional and general graphics functions Two dimensional graphics functions plot Plot lines and points in 2D space fill Draw filled 2D polygons eldraw2 Draw undeformed finite element mesh eldisp2 Draw deformed finite element mesh eldia2 Draw section force diagram elflux2 Plot flux vectors eliso2 Draw isolines for nodal quantities elprinc2 Plot principal stresses General graphics functions axis Axis scaling and appearance clf Clear current figure figure Create figures grid Grid lines hold Hold current graph print Print graph or save graph to file text Add text to current plot title Titles for 2D and 3D plots xlabel Axis labels for 2D and 3D plots ylabel zlabel 8 1 GRAPHICS axis Purpose Plot axis scaling and appearance Syntax axis xmin xmax ymin ymax axis xmin xmax ymin ymax zmin zmax axis auto axis square axis equal axis off axis on Description axis xmin xmax ymin ymax sets scaling for the x and y axes on the current 2D plot axis xmin xmax ymin ymax zmin zmax sets the scaling for the x y and z axes on the current 3D plot axis auto returns the axis scaling to its default automatic mode where for each plot xmin min x xmax max x etc axis square makes the current axis box square i
133. tangular Turner Clough element plantcs Compute stresses and strains plani4e Compute element matrices 4 node isoparametric element plani4ds Compute stresses and strains plani4f Compute internal element forces plani8e Compute element matrices 8 node isoparametric element plani8s Compute stresses and strains plani8f Compute internal element forces 3D solid functions soli8e Compute element matrices 8 node isoparametric element soli8s Compute stresses and strains soli8f Compute internal element forces 5 9 3 ELEMENT plante Two dimensional solid elements Purpose Compute element matrices for a triangular element in plane strain or plane stress Ug U3 Syntax Ke plante ex ey ep D Ke fe plante ex ey ep D eq Description plante provides an element stiffness matrix Ke and an element load vector fe for a triangular element in plane strain or plane stress The element nodal coordinates x1 Y1 2 etc are supplied to the function by ex and ey The type of analysis ptype and the element thickness t are supplied by ep ptype 1 plane stress ptype 2 plane strain and the material properties are supplied by the constitutive matrix D Any arbitrary D matrix with dimensions from 3 x 3 to 6 x 6 may be given For an isotropic elastic material the constitutive matrix can be formed by the function hooke see Section 4 ex 21 T2 z3 ep ptype t eye ae ae p ptype t D D D B 11 12 13 B D31 D32 D
134. te complex frequency response functions Syntax Y sweep K C M p w Description sweep computes the complex frequency response function for a system of the form K iwC w M y w p Here K C and M represent the m by m stiffness damping and mass matrices re spectively The vector p defines the amplitude of the force The frequency response function is computed for the values of w given by the vector w The complex frequency response function is stored in the matrix Y with dimension m by n where n is equal to the number of circular frequencies defined in w Example The steady state response can be computed by gt gt X real Y exp i w t and the amplitude by gt gt Z abs Y SYSTEM 6 3 14 7 Statements and macros Statements describe algorithmic actions that can be executed There are two different types of control statements conditional and repetitive The first group defines conditional jumps whereas the latter one defines repetition until a conditional statement is fulfilled Macros are used to define new functions to the MATLAB or CALFEM structure or to store a sequence of statements in an m file Control statements if Conditional jump for Initiate a loop while Define a conditional loop Macros function Define a new function script Store a sequence of statements STATEMENTS if Purpose Conditional jump Syntax if logical expression els if logical
135. tem functions eigen Purpose Solve the generalized eigenvalue problem Syntax L eigen K M L eigen K M b L X eigen K M L X eigen K M b Description eigen solves the eigenvalue problem K AM 0 where K and M are square matrices The eigenvalues are stored in the vector L and the corresponding eigenvectors in the matrix X If certain rows and columns in matrices K and M are to be eliminated in computing the eigenvalues b must be given in the function The rows and columns that are to be eliminated are described in the vector b defined as dof do fa do fnb The computed eigenvalues are given in order ranging from the smallest to the largest The eigenvectors are normalized in order that XTMX I where I is the identity matrix 6 2 5 SYSTEM extract Static system functions Purpose Extract element nodal quantities from a global solution vector amp 8 u edof eln i j m n i u3 ed u u uz uz a 8 Syntax ed extract edof a Description extract extracts element displacements or corresponding quantities a from the global solution vector a stored in a Input variables are the element topology matrix edof defined in assem and the global solution vector a The output variable ed a contains the element displacement vector If Edof contains more than one element Ed will be a matrix nel where row i gives the element displacements for the e
136. the Jacobian matrix r Ox ae ot On OC ya WY Oy y OE On o az O2 ds ot On OC Evaluation of the integrals is done by Gauss integration 5 4 15 ELEMENT flw3i8s Three dimensional heat flow elements Purpose Compute heat flux and temperature gradients in an 8 node isoparametric heat flow element Syntax es et eci flw3i8s ex ey ez ep D ed Description flw3i8s computes the heat flux vector es and the temperature gradient et or corre sponding quantities in an 8 node isoparametric heat flow element The input variables ex ey ez ep and the matrix D are defined in flw3i8e The vector ed contains the nodal temperatures a of the element and is obtained by the function extract as ed a T To T3 rE T The output variables dz dy qd 2 2 2 avg a E n3 dt YW amp oT aT T Ox Oy Oz OF T rT m1 MW a aa orr oe Oy aa ee OT ar ar ee Ox Oy Oz contain the heat flux the temperature gradient and the coordinates of the integra tion points The index n denotes the number of integration points used within the element cf flw3i8e Theory The temperature gradient and the heat flux are computed according to VT B a q DVT where the matrices D B and a are described in flw3i8e and where the integration points are chosen as evaluation points ELEMENT 5 4 16 5 5 Solid elements Solid elements are available for two dimensional analysi
137. thod Only the cases ptype 2 3 and 4 are implemented 4 3 MATERIAL dmises Purpose Form the elasto plastic continuum matrix for an isotropic hardening von Mises ma terial Syntax D dmises ptype mp es st Description dmises forms the elasto plastic continuum matrix for an isotropic hardening von Mises material The input variable ptype is used to define the type of analysis cf hooke The vector mp contains the material constants mp Evh where F is the modulus of elasticity v is the Poisson s ratio and h is the plastic modulus The matrix es contains current stresses obtained from plants or some similar s function and the vector st contains the current state parameters st yi oy erg where yi 1 if the material behaviour is elasto plastic and yi 0 if the material behaviour is elastic The current yield stress is denoted by o and the current effective plastic strain by epp Note Only the case ptype 2 is implemented MATERIAL 4 4 5 Element functions 5 1 Introduction The group of element functions contains functions for computation of element matrices and element forces for different element types The element functions have been divided into the following groups Spring element Bar elements Heat flow elements Solid elements Beam elements Plate element For each element type there is a function for computation of the element stiffness
138. tress ptype 1 the D matrix is reduced to a 3 x 3 matrix by static condensation using oz Grz Oyz O These stress components are connected with the rows 3 5 and 6 in the D matrix respectively If a larger D matrix than 3 x 3 is used for plane strain ptype 2 the D matrix is reduced to a 3 x 3 matrix using Yez Fiz 0 This implies that a 3 x 3 D matrix is created by the rows and the columns 1 2 and 4 from the original D matrix Evaluation of the integrals is done by Gauss integration ELEMENT 5 5 28 Two dimensional solid elements plani8s Purpose Compute stresses and strains in an 8 node isoparametric element in plane strain or plane stress Syntax es et eci plani8s ex ey ep D ed Description plani8s computes stresses es and the strains et in an 8 node isoparametric element in plane strain or plane stress The input variables ex ey ep and the matrix D are defined in plani8e The vector ed contains the nodal displacements a of the element and is obtained by the function extract as ed a u UQ sas Ure The output variables 1 1 1 1 1 1 Org O jy loz Ony Tz yz 2 2 2 2 2 2 es ot Ore Oyy oz Oxy es Oyz n2 n2 n n n2 n2 Orr Oyy lor Ory Ozz Oyz 1 1 1 1 1 n1 amp E La yy zz Vry Yarz Ty Ly YW 2 2 2 2 2 2 E E T Y2 et eT r yy A Vry Yaz Vyz eci l l n n n n
139. ty EF and the cross section area A Theory The global element stiffness matrix K is computed according to K G K G where EA 1 l N z Nyz N O0 0 0 The transformation matrix G contains the direction cosines _ ty T m Y2 y Z2 1 T L yT L z L where the length L y 22 1 yo y1 z2 21 ELEMENT 5 3 6 Three dimensional bar element bar3s Purpose Compute normal force in a three dimensional bar element Syntax es bar3s ex ey ez ep ed Description bar3s computes the normal force in a three dimensional bar element The input variables ex ey ez and ep are defined in bar3e and the element nodal displacements stored in ed are obtained by the function extract The output variable es N contains the normal force N of the bar Theory The normal force N is computed from EA N T 1 1 G af where E A L and the transformation matrix G are defined in bar3e The nodal displacements in global coordinates e T af u U2 ug u4 Us usg are also shown in bar3e Note that the transpose of a is stored in ed 5 3 7 ELEMENT 5 4 Heat flow elements Heat flow elements are available for one two and three dimensional analysis For one dimensional heat flow the spring element spring1 is used A variety of important physical phenomena are described by the same differential equa tion as the heat flow problem The heat flow element is thus
140. ust be included The output variable xf uM No Ve Ms contains the section forces at the ends of the beam ELEMENT 5 6 18 Two dimensional beam element beam2gs Theory The section forces at the ends of the beam are obtained from the element force vector P N V M N V M computed according to P K Ga f The matrix G is described in beam2e The matrix K and the nodal displacements e T af u U2 u3 Us Us use are described in beam2g Note that the transpose of af is stored in ed 5 6 19 ELEMENT beam2d Two dimensional beam element Purpose Compute element stiffness mass and damping matrices for a two dimensional beam element XY y ee R Uy xpy Syntax Ke Me beam2d ex ey ep Ke Me Ce beam2d ex ey ep Description beam2d provides the global element stiffness matrix Ke the global element mass matrix Me and the global element damping matrix Ce for a two dimensional beam element The input variables ex and ey are described in beam2e and ep E AI m a a l contains the modulus of elasticity the cross section area A the moment of inertia I the mass per unit length m and the Raleigh damping coefficients ap and a1 If ag and a are omitted the element damping matrix Ce is not computed ELEMENT 5 6 20 Two dimensional beam element beam2d Theory The element stiffness matrix K the element mass matrix M and the element damping
141. utput vector in this case will be equal to N See also the description of the command fft The inverse Fourier coefficients x j stored in the variable x are computed according to ulk VE Mz z j 1 N k 1 where wy eo emt N Note This is a MATLAB built in function See also fft 6 3 7 SYSTEM ritz Dynamic system functions Purpose Compute approximative eigenvalues and eigenvectors by the Lanczos method Syntax L ritz K M f m L ritz K M f m b L X ritz K M f m L X ritz K M f m b Description ritz computes by the use of the Lanczos algorithm m approximative eigenvalues and m corresponding eigenvectors for a given pair of n by n matrices K and M and a given non zero starting vector f If certain rows and columns in matrices K and M are to be eliminated in computing the eigenvalues b must be given in the command The rows and columns to be eliminated are described in the vector b defined as dof do fo do fnb Note If the number of vectors m is chosen less than the total number of degrees of freedom n only about the first m 2 Ritz vectors are good approximations of the true eigenvectors Recall that the Ritz vectors satisfy the M orthonormality condition X MX I where I is the identity matrix SYSTEM 6 3 8 Dynamic system functions spectra Purpose Compute seismic response spectra for elastic design Syntax s spectra a xi dt f Description spect
142. vector x The commands title text xlabel zlabel ylabel ylabel write text as a title of the current plot and label and ylabel as labels of the coordinate axis Grid lines are added with grid and clf clears the current figure EXAMPLES 9 2 12 9 3 Static analysis This section illustrates some linear static finite element calculations The examples deal with structural problems as well as field problems such as heat conduction Static analysis exs1 exs2 exs3 exs4 exs5 exs6 exs Linear spring system One dimensional heat flow Simply supported beam Plane truss Plane frame Geometry based frame analysis Two dimensional diffusion The introductory example exs1 illustrates the basic steps in the finite element method for a simple structure of one dimensional linear springs The linear spring or analogy element is also used in example exs2 to solve a problem of heat conduction through a wall A simply supported beam is analysed in example exs3 Element forces and reactions at the supports are calculated A two dimensional plane truss is analysed in example exs4 and a two dimensional plane frame is analysed in example exs5 A frame built up from both beams and bars is analysed in example exs6 The finite element model is defined from the geometry Graphics facilities are also explained in example exs5 and exs6 Note The examples listed above are supplied as m files on the
143. ven inside the parentheses and also separated by commas The statement help function provides information about purpose and syntax for the specified function 1 1 INTRODUCTION The available functions are organized in groups as follows Each group is described in a separate chapter Groups of functions General purpose commands for managing variables workspace output etc Matrix functions for matrix handling Material functions for computing material matrices Element functions for computing element matrices and element forces System functions for setting up and solving systems of equations Statement functions for algorithm definitions Graphics functions for plotting INTRODUCTION 1 2 2 General purpose functions The general purpose functions are used for managing variables and workspace control of output etc The functions listed here are a subset of the general purpose functions described in the MATLAB manual The functions can be divided into the following groups Managing commands and functions help Online documentation type List m file what Directory listing of m mat and mex files ue Continuation Write a comment line Managing variables and the workspace clear Remove variables from workspace disp Display variables in workspace on display screen load Retrieve variable from disk and load in workspace save Save matrix bank variable on disk who List directory of varia
144. www byggmek 1th se Calfem Preface CALFEM is an interactive computer program for teaching the finite element method FEM The name CALFEM is an abbreviation of Computer Aided Learning of the Finite Element Method The program can be used for different types of structural mechanics problems and field problems CALFEM the program and its built in philosophy have been developed at the Division of Structural Mechanics starting in the late 70 s Many coworkers former and present have been engaged in the development at different stages of whom we might mention Per Erik Austrell Hakan Carlsson Ola Dahlblom Jonas Lindemann Anders Olsson Karl Gunnar Olsson Kent Persson Anders Peterson Hans Petersson Matti Ristinmaa Goran Sandberg The present release of CALFEM as a toolbox to MATLAB represents the latest develop ment of CALFEM The functions for finite element applications are all MATLAB functions M files as described in the MATLAB manual We believe that this environment increases the versatility and handling of the program and above all the ease of teaching the finite element method Lund November 22 2000 Division of Structural Mechanics and Division of Solid Mechanics Contents 1 Introduction 1 1 2 General purpose functions 2 1 3 Matrix functions 3 1 4 Material functions 4 1 5 Element functions 5 1 1 pd Antroduction z cs 4k n aA ea A e AA a oe te A 5 1 1 92 OPNE lement a ia tra R a aoe eR Gf i aae EAE ie iaia
145. x Edof a global stiffness matrix K and load vector f are defined The element matrices Ke and fe are computed by the function beam2e These matrices are then assembled in the global matrices using the function assem gt gt Edof 1 123456 2 101112789 3 456789 gt gt K zeros 12 f zeros 12 1 4 1000 gt gt A1l 45 3e 4 A2 142 8e 4 EXAMPLES 9 3 16 Static analysis exs5 gt gt 1 1 2510e 8 I12 33090e 8 gt gt E 2 1e11 gt gt ep1 E A1 I1 ep3 E A2 I2 gt gt exi 0 0 ex2 6 6 ex3 0 6 gt gt ey1 0 4 ey2 0 4 ey3 4 4 gt gt eqi 0 0 gt gt eq2 0 0 gt gt eq3 0 75000 gt gt Kel beam2e ex1 eyl ep1 gt gt Ke2 beam2e ex2 ey2 ep1 gt gt Ke3 fe3 beam2e ex3 ey3 ep3 eq3 gt gt K assem Edof 1 K Ke1 gt gt K assem Edof 2 K Ke2 gt gt K f assem Edof 3 K Ke3 f fe3 The system of equations are solved considering the boundary conditions in be gt gt bc 1 0 2 0 3 0 10 0 11 0 12 0 gt gt a solveq K f bc 0 0006 0 0009 0 0079 0 0005 0 0009 0 0079 9 3 17 EXAMPLES exs5 Static analysis The element displacements are obtained from the function extract and the function beam2s computes the section forces gt gt Ed extract Edof a gt gt esl beam2s ex1 ey1 ep1 Fd 1 eq1 20 esi 1 0e 05 2 2467 0 1513 0 1981 2 2467 0 1513 0 1662 2 2467 0 1513 0 4070 gt gt es2 beam2s ex2 ey2
146. x1 Fy1 1 3 1 eldraw2 Ex2 Fy2 1 2 1 The element stiffness matrices are generated and assembled in two loops one for the beams and one for the bars The element stiffness functions beam2e and bar2e use the element coordinate matrices ex and ey These matrices are extracted from the global coordinates Coord by the function coordxtr above EXAMPLES 9 3 24 Static analysis exs6 for i 1 6 Ke beam2e Ex1 i Ey1 i ep1 K assem Edofi i K Ke end for i 1 4 Ke bar2e Ex2 i Ey2 i ep2 K assem Edof2 i K Ke end The global system of equations is solved considering the boundary conditions in bc bc 1 0 20 3 0 40 5 0 60 a Q solveq K f bc and the deformed frame is displayed by the function eldisp2 where the displacements are scaled by the variable magnfac Edi extract Edof1 a Ed2 extract Edof2 a magnfac eldisp2 Ex1 Ey1 Ed1 eldisp2 Ex2 Fy2 Ed2 2 1 1 magnfac grid 9 3 25 EXAMPLES Static analysis exs6 2H T re sess ey Taa D l I5 P ke 1 l 7 l 4 l v I No Z l SK I NG r I ER l EEIN E SENE Sten 1 5 I E f I 7 I man I A I Z be of Ne va ce WZ lt i ENNET TEETE zos q 1 E lt gt Z jl ni l X A I I 7 I 4 A NVA OS sets Pesce aatineah EANET py creas CEEE SE J N Za X A X do X 0 L 1 0 5 0 0 5 1 1 5 x EXAMPLES 9 3 26 Static analysis exs7 Purpose Analysis of tw
147. x2 Ex Ey Es plotpar magnfac Description elflux2 displays element heat flux vectors or corresponding quantities for a number of elements of the same type The flux vectors are displayed as arrows at the element centroids Note that only the flux vectors are displayed To display the element mesh use eldraw2 Input variables are the coordinate matrices Ex and Ey and the element flux matrix Es defined in flw2ts or flw2qs The variable plotpar sets plot parameters for the flux arrows plotpar arrowtype arrowcolor arrowtype 1 solid arrowcolor 1 black 2 dashed 2 blue 3 dotted 3 magenta 4 red Default if plotpar is omitted is solid black arrows The magnification factor magnfac is a scalar that magnifies the arrows in relation to the element size The magnification factor is set automatically if it is omitted in the input list Limitations Supported elements are triangular 3 node and quadrilateral 4 node elements GRAPHICS 8 8 eliso2 Purpose Display element iso lines for two dimensional elements Syntax eliso2 Ex Ey Ed isov eliso2 Ex Ey Ed isov plotpar Description eliso2 displays element iso lines for a number of elements of the same type Note that only the iso lines are displayed To display the element mesh use eldraw2 Input variables are the coordinate matrices Ex and Ey formed by the function co ordxtr and the element nodal quantities e g displacement or energy potential

Calfem - Chair of Applied Mathematics / Numerical Analysis

Contents

Download Pdf Manuals

Related Search

Related Contents