Contents
1 Functio unctions ns — By cate categor gory y
1
1.1 Genera Generall functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2 Postp Postproce rocessi ssing ng . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.3 1.3
Dyna Dynami mics cs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
1.4 1.4
Beam Be am funct functio ions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
1.5 Truss functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
1.6 Genera Generall shel shelll functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
1.7 Shell4 Shell4 functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
1.8 Shell8 Shell8 functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
1.9 Nonline Nonlinear ar anal analysi ysiss functi functions ons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
2 Tutor utoria iall
71
2.1 Tutorial utorial:: static static anal analysi ysiss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
2.2 Tutorial utorial:: dynamic dynamic analysi analysiss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
2.2. 2.2.11
Eige Eigenm nmode odess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
2.2.2
Modal superposition superposition:: time domain: domain: piecewise piecewise exact exact integrat integration ion . . . . . . . . .
79
2.2.3 2.2 .3
Modal superpos superpositi ition: on: transform transform to freq frequen uency cy domai domain n. . . . . . . . . . . . . . .
81
2.2.4 2.2 .4
Direct Direct time time integ integrat ration ion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
2.2.5 2.2 .5
Direct Direct solu solutio tion n in the frequenc frequency y domain domain . . . . . . . . . . . . . . . . . . . . . .
85
2.2. 2.2.66
Compa Compari riso son n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
86
2.3 Tutorial utorial:: shell analys analysis is . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
88
2.4 Tutorial utorial:: nonline nonlinear ar anal analysi ysiss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
93
3 Ex Exam ampl ples es
97
3.1 3.1
Exam Exampl plee 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
98
3.2 3.2
Exam Exampl plee 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
3.3 3.3
Exam Exampl plee 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.4 3.4
Exam Exampl plee 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Contents
1 Functio unctions ns — By cate categor gory y
1
1.1 Genera Generall functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2 Postp Postproce rocessi ssing ng . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.3 1.3
Dyna Dynami mics cs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
1.4 1.4
Beam Be am funct functio ions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
1.5 Truss functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
1.6 Genera Generall shel shelll functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
1.7 Shell4 Shell4 functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
1.8 Shell8 Shell8 functio functions ns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
1.9 Nonline Nonlinear ar anal analysi ysiss functi functions ons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
2 Tutor utoria iall
71
2.1 Tutorial utorial:: static static anal analysi ysiss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
2.2 Tutorial utorial:: dynamic dynamic analysi analysiss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
2.2. 2.2.11
Eige Eigenm nmode odess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
2.2.2
Modal superposition superposition:: time domain: domain: piecewise piecewise exact exact integrat integration ion . . . . . . . . .
79
2.2.3 2.2 .3
Modal superpos superpositi ition: on: transform transform to freq frequen uency cy domai domain n. . . . . . . . . . . . . . .
81
2.2.4 2.2 .4
Direct Direct time time integ integrat ration ion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
2.2.5 2.2 .5
Direct Direct solu solutio tion n in the frequenc frequency y domain domain . . . . . . . . . . . . . . . . . . . . . .
85
2.2. 2.2.66
Compa Compari riso son n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
86
2.3 Tutorial utorial:: shell analys analysis is . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
88
2.4 Tutorial utorial:: nonline nonlinear ar anal analysi ysiss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
93
3 Ex Exam ampl ples es
97
3.1 3.1
Exam Exampl plee 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
98
3.2 3.2
Exam Exampl plee 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
3.3 3.3
Exam Exampl plee 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.4 3.4
Exam Exampl plee 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
ii
Contents
3.5 3.5
Exam Exampl plee 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
3.6 3.6
Exam Exampl plee 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
3.7 3.7
Exam Exampl plee 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
4 Functions unctions — Alpha Alphabetica beticall list list
127
Chapter 1
Functions — By category 1.1
General functions
getdof asmkm removedof addconstr tconstr nodalvalues elemloads accel elemforces elemdisp selectdof unselectdof selectnode reprow multdloads
Get the vector with the degrees of freedom of the model. Assemble stiffness and mass matrix. Remove DOF with Dirichlet boundary conditions equal to zero. Add constraint equations to the stiffness matrix and load vector. Return matrices to apply constraint equations. Construct a vector with the values at the selected DOF. Equivalent nodal forces. Compute the distributed loads due to an acceleration. Compute the element forces. Select the element displacements from the global displacement vector. Select degrees of freedom. Unselect degrees of freedom. Select nodes by location. Replicate rows from a matrix. Combine distributed loads.
p.2 p.2 p.2 p.3 p.3 p.4 p.4 p.4 p.5 p.5 p.6 p.7 p.7 p.7 p.8
2
Functions — By category
getdof GETDOF
Get the vector with the degrees of freedom of the model.
DOF=getdof(Elements,Types) builds the vector with the labels of the degrees of freedom for which stiffness is present in the finite element model. Elements Types DOF
Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Degrees of freedom (nDOF * 1)
See also DOF_TRUSS, DOF_BEAM, GETDOF.
removedof REMOVEDOF
Remove DOF with Dirichlet boundary conditions equal to zero.
DOF=removedof(DOF,seldof) removes the specified Dirichlet boundary conditions from the degrees of freedom vector. DOF seldof
Degrees of freedom (nDOF * 1) Dirichlet boundary conditions equal to zero
[NodID.dof]
See also GETDOF.
asmkm ASMKM
Assemble stiffness and mass matrix.
[K,M] = asmkm(Nodes,Elements,Types,Sections,Materials,DOF) K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF) K = asmkm(Nodes,Elements,Types,Sections,Materials) assembles the stiffness and the mass matrix using the finite element method. Nodes Elements Types Sections Materials DOF K
Node definitions [NodID Element definitions [EltID Element type definitions {TypID Section definitions [SecID Material definitions [MatID Degrees of freedom (nDOF * 1) Stiffness matrix (nDOF * nDOF)
x y z] TypID SecID MatID n1 n2 ...] EltName Option1 ... } SecProp1 SecProp2 ...] MatProp1 MatProp2 ... ]
General functions
M
Mass matrix (nDOF * nDOF)
See also KE_TRUSS, KE_BEAM.
addconstr ADDCONSTR
Add constraint equations to the stiffness matrix and load vector.
[K,F]=addconstr(Constr,DOF,K,F) [K,F,M]=addconstr(Constr,DOF,K,[],M) [K,F,M]=addconstr(Constr,DOF,K,F,M) modifies the stiffness matrix, the mass matrix and the load vector according to the applied constraint equations. The dimensions of the stiffness matrix, the mass matrix and the load vector are kept the same. The resulting stiffness and mass matrix are not symmetric anymore. This function can be used as well to apply imposed displacements. Constr
DOF K F M
Constraint equation: Constant=CoefS*SlaveDOF+CoefM1*MasterDOF1+CoefM2*MasterDOF2+... [Constant CoefS SlaveDOF CoefM1 MasterDOF1 CoefM2 MasterDOF2 ...] Degrees of freedom (nDOF * 1) Stiffness matrix (nDOF * nDOF) Load vector (nDOF * nSteps) Mass matrix (nDOF * nDOF)
tconstr TCONSTR
Return matrices to apply constraint equations.
[T,Q0,MasterDOF]=tconstr(Constr,DOF) returns matrices to apply constraint equations to the stiffness and mass matrix and the load vector: Kr=T.’*K*T, Mr=T.’*M*T and Fr=T.’*(F-K*Q0). The original displacement vector is computed using U=T*Ur+Q0. Constr
DOF
Constraint equation: Constant=CoefS*SlaveDOF+CoefM1*MasterDOF1+CoefM2*MasterDOF2+... [Constant CoefS SlaveDOF CoefM1 MasterDOF1 CoefM2 MasterDOF2 ...] Degrees of freedom (nDOF * 1)
3
4
Functions — By category
nodalvalues NODALVALUES
Construct a vector with the values at the selected DOF.
F=nodalvalues(DOF,seldof,values) constructs a vector with the values at the selected DOF. This function can be used to obtain a load vector, initial displacements, velocities or accelerations. DOF seldof values F
Degrees of freedom (nDOF * 1) Selected degrees of freedom [NodID.dof] (nValues * 1) Corresponding values [Value] (nValues * nSteps) Load vector (nDOF * nSteps)
See also ELEMLOADS.
elemloads ELEMLOADS
Equivalent nodal forces.
F=elemloads(DLoads,Nodes,Elements,Types,DOF) computes the equivalent nodal forces of a distributed load (in the global coordinate system). DLoads Nodes Elements Types DOF F
Distributed loads [EltID n1globalX n1globalY n1globalZ ...] Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Degrees of freedom (nDOF * 1) Load vector (nDOF * 1)
See also LOADS_TRUSS, LOADS_BEAM, NODALVALUES.
accel ACCEL
Compute the distributed loads due to an acceleration.
DLoads=accel(Accelxyz,Elements,Types,Sections,Materials) computes the distributed loads due to an acceleration. In order to simulate gravity, accelerate the structure in the direction opposite to gravity.
General functions
Accelxyz Elements Types Sections Materials DLoads
Acceleration [Ax Ay Az] (1 * 3) Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Section definitions [SecID SecProp1 SecProp2 ...] Material definitions [MatID MatProp1 MatProp2 ... ] Distributed loads [EltID n1globalX n1globalY n1globalZ ...]
See also ELEMLOADS, ACCEL_BEAM, ACCEL_TRUSS.
elemforces ELEMFORCES
Compute the element forces.
[ForcesLCS,ForcesGCS] =elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads) [ForcesLCS,ForcesGCS] =elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U) computes the element forces in the local (beam convention) and the global (algebraic convention) coordinate system. Nodes Elements Types Sections Materials DOF U DLoads ForcesLCS
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Section definitions [SecID SecProp1 SecProp2 ...] Material definitions [MatID MatProp1 MatProp2 ... ] Degrees of freedom (nDOF * 1) Displacements (nDOF * 1) Distributed loads [EltID n1globalX n1globalY n1globalZ ...] Element forces in LCS (beam convention) [N Vy Vz T My Mz] (nElem * 12) ForcesGCS Element forces in GCS (algebraic convention) (nElem * 12) See also FORCES_TRUSS, FORCES_BEAM.
elemdisp ELEMDISP
Select the element displacements from the global displacement vector.
UeGCS = elemdisp(Type,NodeNum,DOF,U) selects the element displacements from the global displacement vector. Type
Element type e.g. ’beam’,’truss’, ...
5
6
Functions — By category
M yj
x z
T j
j M zj
M yi
V zj
V zi
T i
N i
N j
V yj
y
i
M zi
V yi
Figure 1.1: Beam convention. NodeNum DOF U UeGCS
Node numbers (1 * nNodes) Degrees of freedom (nDOF * 1) Displacements (nDOF * 1) Elements displacements
See also ELEMFORCES, DOF_BEAM, DOF_TRUSS.
selectdof SELECTDOF
Select degrees of freedom.
L=selectdof(DOF,seldof) [L,I]=selectdof(DOF,seldof) L=selectdof(DOF,seldof,’Ordering’,ordering) creates the matrix to extract degrees of freedom from the global degrees of freedom by matrix multiplication. DOF seldof L I ordering
Degrees of freedom (nDOF * 1) Selected DOF labels (kDOF * 1) Selection matrix (kDOF * nDOF) Index vector (kDOF * 1) ’seldof’,’DOF’ or ’sorted’
Default: ’seldof’
General functions
Ordering of L and I similar as seldof, DOF or sorted
unselectdof UNSELECTDOF
Unselect degrees of freedom.
L=unselectdof(DOF,seldof) [L,I]=unselectdof(DOF,seldof) creates the matrix to unselect degrees of freedom from the global degrees of freedom. DOF seldof L I
Degrees of freedom (nDOF * 1) Unselected dof labels (ndof * 1) Selection matrix ((nDOF-ndof) * nDOF) Index vector ((nDOF-ndof) * 1)
See also SELECTDOF.
selectnode SELECTNODE
Select nodes by location.
Nodesel=selectnode(Nodes,x,y,z) Nodesel=selectnode(Nodes,xmin,ymin,zmin,xmax,ymax,zmax) selects nodes by location. Nodes Nodesel
Node definitions [NodID x y z] Node definitions of the selected nodes
reprow REPROW
Replicate rows from a matrix.
Matrix=reprow(Matrix,RowSel,nTime,RowInc) replicates the selected rows from a matrix a number of times and adds them below the existing rows. The k-th time the increment RowInc is added k times to the copied rows.
7
8
Functions — By category
Matrix RowSel nTime RowInc
Matrix (nRow * nCol) Rows to be copied (1 * kRow) Number of times Increments that are added to the different columns of the copied rows (1 * nCol)
multdloads MULTDLOADS
Combine distributed loads.
DLoads=multdloads(DLoads_1,DLoads_2,...,DLoads_k) combines the distributed loads of multiple load cases into one 3D array. Each plane corresponds to one load case. DLoads_k
Distributed loads
DLoads
Distributed loads
See also ELEMLOADS.
[EltID n1globalX n1globalY n1globalZ ...] (nElem_k * 7) [EltID n1globalX n1globalY n1globalZ ...] (maxnElem * 7 * k)
Postprocessing
1.2
9
Postprocessing
plotnodes Plot the nodes. Plot the elements. plotelem plotdisp Plot the displacements. plotforc Plot the forces. plotlcs Plot the local element coordinate systems. plotstress Plot the stresses. Plot filled contours of stresses. plotstresscontourf animdisp Animate the displacements. getmovie Get the movie from a figure where an animation has been played. printdisp Display the displacements in the command window. printforc Display the forces in the command window.
p.10 p.10 p.11 p.11 p.12 p.13 p.14 p.14 p.15 p.15 p.16
10
Functions — By category
plotnodes PLOTNODES
Plot the nodes.
plotnodes(Nodes) plots the nodes. Nodes
Node definitions
[NodID x y z]
plotnodes(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’Numbering’ Plots the node numbers. Default: ’on’. ’GCS’ Plots the global coordinate system. Default: ’on’. ’Handle’ Plots in the axis with this handle. Default: current axis. Additional parameters are redirected to the PLOT3 function which plots the nodes. h = PLOTNODES(...) returns a struct h with handles to all the objects in the plot. See also PLOTELEM, PLOTDISP.
plotelem PLOTELEM
Plot the elements.
plotelem(Nodes,Elements,Types) plots the elements. Nodes Elements Types
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... }
plotelem(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’Numbering’ Plots the element numbers. Default: ’on’. ’GCS’ Plots the global coordinate system. Default: ’on’. ’Handle’ Plots in the axis with this handle. Default: current axis. Additional parameters are redirected to the PLOT3 function which plots the elements. h = PLOTELEM(...) returns a struct h with handles to all the objects in the plot. See also COORD_TRUSS, COORD_BEAM, PLOTDISP.
Postprocessing
11
plotdisp PLOTDISP
Plot the displacements.
plotdisp(Nodes,Elements,Types,DOF,U,DLoads,Sections,Materials) plotdisp(Nodes,Elements,Types,DOF,U,[],Sections,Materials) plotdisp(Nodes,Elements,Types,DOF,U) DispScal=plotdisp(Nodes,Elements,Types,DOF,U,DLoads,Sections,Materials) plots the displacements. If DLoads, Sections and Materials are supplied, the displacements that occur due to the distributed loads if all nodes are fixed, are superimposed. Nodes Elements Types DOF U DLoads
Sections Materials DispScal
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Degrees of freedom (nDOF * 1) Displacements (nDOF * 1) Distributed loads [EltID n1globalX n1globalY n1globalZ ...] (use empty array [] when shear deformation (in beam element) is considered but no DLoads are present) Section definitions [SecID SecProp1 SecProp2 ...] Material definitions [MatID MatProp1 MatProp2 ... ] Displacement scaling
plotdisp(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’DispScal’ Displacement scaling. Default: ’auto’. ’DispMax’ Mention maximal displacement. Default: ’on’. ’Undeformed’ Plots the undeformed mesh. Default: ’k:’. ’Handle’ Plots in the axis with this handle. Default: current axis. Additional parameters are redirected to the PLOT3 function which plots the deformations. [DispScal,h] = plotdisp(...) returns a struct h with handles to all the objects in the plot. See also DISP_TRUSS, DISP_BEAM, PLOTELEM.
plotforc PLOTFORC
Plot the forces.
plotforc(ftype,Nodes,Elements,Types,Forces,DLoads) plotforc(ftype,Nodes,Elements,Types,Forces) plots the forces (in beam convention).
12
Functions — By category
ftype
Nodes Elements Types Forces DLoads
’norm’ Normal force (in the local x-direction) ’sheary’ Shear force in the local y-direction ’shearz’ Shear force in the local z-direction ’momx’ Torsional moment (around the local x-direction) ’momy’ Bending moment around the local y-direction ’momz’ Bending moment around the local z-direction Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Element forces in LCS (beam convention) [N Vy Vz T My Mz] (nElem * 12) Distributed loads [EltID n1globalX n1globalY n1globalZ ...]
plotforc(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’ForcScal’ Force scaling. Default: ’auto’. ’Values’ Force values. Default: ’on’. ’Undeformed’ Plots the undeformed mesh. Default: ’k-’. ’Handle’ Plots in the axis with this handle. Default: current axis. Additional parameters are redirected to the PLOT3 function which plots the forces. [ForcScal,h] = plotforc(...) returns a struct h with handles to all the objects in the plot. See also FDIAGRGCS_BEAM, FDIAGRGCS_TRUSS.
plotlcs PLOTLCS
Plot the local element coordinate systems.
[h,vLCS] = plotlcs(Nodes,Elements,Types) [h,vLCS] = plotlcs(Nodes,Elements,Types,[],varargin) plotlcs(Nodes,Elements,Types,vLCS,varargin) Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } vLCS Element coordinate systems (nElem * 9) plotlcs(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’GCS’ Plot the GCS. Default: ’on’. ’Undeformed’ Plots the undeformed mesh. Default: ’k-’. ’Handle’ Plots in the axis with this handle. Default: current axis.
Postprocessing
13
See also ELEMSTRESS
plotstress PLOTSTRESS
Plot the stresses.
plotstress(stype,Nodes,Elements,Types,Sections,Forces,DLoads) plotstress(stype,Nodes,Elements,Types,Sections,Forces) plots the stresses. stype
’snorm’ ’smomyt’
Nodes Elements Types Sections Forces
Normal stress due to normal force Normal stress due to bending moment around the local y-direction at the ’smomyb’ Normal stress due to bending moment around the local y-direction at the ’smomzt’ Normal stress due to bending moment around the local z-direction at the ’smomzb’ Normal stress due to bending moment around the local z-direction at the ’smax’ Maximal normal stress (normal force ’smin’ Minimal normal stress (normal force Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID Element type definitions {TypID EltName Option1 Section definitions [A ky kz Ixx Iyy Izz yt Element forces in LCS (beam convention) [N Vy Vz
DLoads
Distributed loads
top bottom top bottom and bending moment) and bending moment)
n1 n2 ...] ... } yb zt zb] T My Mz] (nElem * 12) [EltID n1globalX n1globalY n1globalZ ...]
plotstress(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’StressScal’ Stress scaling. Default: ’auto’. ’Values’ Stress values. Default: ’on’. ’Undeformed’ Plots the undeformed mesh. Default: ’k-’. ’Handle’ Plots in the axis with this handle. Default: current axis. Additional parameters are redirected to the PLOT3 function which plots the stresses. [StressScal,h] = plotstress(...) returns a struct h with handles to all the objects in the plot. See also SDIAGRGCS_BEAM, SDIAGRGCS_TRUSS.
14
Functions — By category
plotstresscontourf PLOTSTRESSCONTOURF
Plot filled contours of stresses.
plotstresscontourf(stype,Nodes,Elements,Types,S) plotstresscontourf(stype,Nodes,Elements,Types,S,DOF,U) plots stress contours (output from ELEMSTRESS). stype
’sx’ Normal stress (in the global/local x-direction) ’sy’ ’sz’ ’sxy’ Shear stress ’syz’ ’sxz’ Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } S Element stresses in GCS/LCS [sxx syy szz sxy syz sxz] (nElem * 72) plotstresscontourf(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’location’ Location (top,mid,bot). Default: ’top’. ’GCS’ Plot the GCS. Default: ’on’. ’minmax’ Add location of min and max stress. Default: ’on’. ’colorbar’ Add colorbar. Default: ’on’. ’ncolor’ Number of colors in colormap. Default: 10. ’Handle’ Plots in the axis with this handle. Default: current axis. See also ELEMSTRESS.
animdisp ANIMDISP
Animate the displacements.
DispScal=animdisp(Nodes,Elements,Types,DOF,U) animates the displacements. Nodes Elements Types DOF U DispScal
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Degrees of freedom (nDOF * 1) Displacements (nDOF * nSteps) Displacement scaling
ANIMDISP(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified:
Postprocessing
15
’DispScal’ ’Handle’ ’Fps’ ’CreateMovie’
Displacement scaling. Default: ’auto’. Plots in the axis with this handle. Default: current axis. Frames per second. Default: 12. Saves the movie in the userdata of the axis of the figure. Use getmovie to get the movie from the axis. Default: ’off’. ’Counter’ Displays the number of the frame for transient displacements. Default: ’on’. Additional parameters are redirected to the PLOTDISP function which plots the individual frames of the movie. See also GETMOVIE, PLOTDISP.
getmovie GETMOVIE
Get the movie from a figure where an animation has been played.
mov=getmovie(h) gets the movie from the userdata of the axis of a figure where an animation has been played. In order to save the movie in the userdata the animation should have been played using animdisp(...,’CreateMovie’,’on’). This function blocks the command prompt until the movie has become available. h mov
Handle of the axis e.g. gca Struct with the frames of the movie
The movie can be played with movie(gcf,mov). See also ANIMDISP, MOVIE.
printdisp PRINTDISP
Display the displacements in the command window.
printdisp(Nodes,DOF,U) displays the displacements in the command window. Nodes DOF U
Node definitions [NodID x y z] Degrees of freedom (nDOF * 1) Displacements (nDOF * 1)
See also PRINTFORC, PLOTDISP.
16
Functions — By category
printforc PRINTFORC
Display the forces in the command window.
printforc(Elements,Forces) displays the forces in the command window. Elements Forces
Element definitions Element forces
See also PRINTDISP.
[EltID TypID SecID MatID n1 n2 ...] [N Vy Vz T My Mz] (nElem * 12)
Dynamics
1.3
17
Dynamics
eigfem msupt msupf cdiff newmark wilson ldlt
Compute the eigenmodes and eigenfrequencies of the finite element model. Modal superposition in the time domain. Modal superposition in the frequency domain. Direct time integration for dynamic systems - central diff. method. Direct time integration for dynamic systems - Newmark method Direct time integration for dynamic systems - Wilson-theta method LDL Gauss factorization of symmetric matrices.
p.18 p.18 p.18 p.19 p.19 p.20 p.20
18
Functions — By category
eigfem EIGFEM
Compute the eigenmodes and eigenfrequencies of the finite element model.
[phi,omega]=eigfem(K,M,nMode) [phi,omega]=eigfem(K,M) computes the eigenmodes and eigenfrequencies of the finite element model. K M nMode phi omega
Stiffness matrix (nDOF * nDOF) Mass matrix (nDOF * nDOF) Number of eigenmodes and eigenfrequencies (default: all) Eigenmodes (in columns) (nDOF * nMode) Eigenfrequencies [rad/s] (nMode * 1)
msupt MSUPT
Modal superposition in the time domain.
x = MSUPT(omega,xi,t,pm,x1,y1,interp) calculates the modal displacements x(t) of a dynamic system with the eigenfrequencies omega and the modal damping ratios xi due to the modal excitation pm(t), given the initial conditions x1 and y1. The solution of the modal differential equations is performed by means of the piecewise exact method. Interp can be ’foh’ (first order hold) or ’zoh’ (zero order hold). In the first case the excitation is assumed to vary linearly within each time step while in the second case it is assumed to be constant: if t(k) <= t < t(k+1) then p(t) = p(k). omega xi t pm x1 y1 interp x
Eigenfrequencies [rad/s] (nMode * 1) Modal damping ratios [ - ] (nMode * 1) or (1 * 1), constant modal damping is assumed in the (1 * 1) case. Time points (1 * N) of the sampling of p, x and t. Modal excitation (nMode * N). Modal displacements at time point t(1) (nMode * 1), defaults to zero. Modal velocities at time point t(1) (nMode * 1), defaults to zero. Interpolation scheme: ’foh’ or ’zoh’, defaults to ’foh’. Modal displacements (nMode * N).
msupf MSUPF
Modal superposition in the frequency domain.
Dynamics
19
[X,H] = MSUPF(omega,xi,Omega,Pm) calculates the modal transfer functions H and the modal displacements x(t) = X * EXP(i * Omega * t) of a dynamic system with the eigenfrequencies omega and the modal damping ratios xi due to the modal excitation pm(t) = Pm * EXP(i * Omega * t). omega xi Omega Pm X H
Eigenfrequencies [rad/s] (nMode * 1) Modal damping ratios [ - ] (nMode * 1) or (1 * 1), constant modal damping is assumed in the (1 * 1) case. Excitation frequencies [rad/s] (1 * N). Complex amplitude of the modal excitation (nMode * N). Complex amplitude of the modal displacements (nMode * N). Modal transfer functions (nMode * N).
cdiff CDIFF Direct time integration for dynamic systems - central diff. method. [u,t] = CDIFF(M,C,K,dt,p,u0,u1) applies the central difference method for the calculation of the nodal displacements u of the dynamic system with the system matrices M, C and K due to the excitation p. M C K dt p u0 u1 u t
Mass matrix (nDof * nDof) Damping matrix (nDof * nDof) Stiffness matrix (nDof * nDof) Time step of the integration scheme (1 * 1). Should be small enough to ensure the stability and the precision of the integration scheme. Excitation (nDof * N). p(:,k) corresponds to time point t(k). Displacements at time point t(1)-dt (nDof * 1). Defaults to zero. Displacements at time point t(1) (nDof * 1). Defaults to zero. Displacements (nDof * N). u(:,k) corresponds to time point t(k). Time axis (1 * N), defined as t = [0:N-1] * dt.
newmark NEWMARK Direct time integration for dynamic systems - Newmark method [u,v,a,t] = NEWMARK(M,C,K,dt,p,u0,v0,a0,[alpha delta]) applies the Newmark method for the calculation of the nodal displacements u, velocities v and accelerations a of the dynamic system with the system matrices M, C and K due to the excitation p. M C K
Mass matrix (nDof * nDof) Damping matrix (nDof * nDof) Stiffness matrix (nDof * nDof)
20
Functions — By category
dt p u0 v0 a0 u t
Time step of the integration scheme (1 * 1). Should be small enough to ensure the stability and the precision of the integration scheme. Excitation (nDof * N). p(:,k) corresponds to time point t(k). Displacements at time point t(1)-dt (nDof * 1). Defaults to zero. Velocities at time point t(1)-dt (nDof * 1). Defaults to zero. Accelerations at time point t(1)-dt (nDof * 1). Defaults to zero. Displacements (nDof * N). u(:,k) corresponds to time point t(k). Time axis (1 * N), defined as t = [0:N-1] * dt.
wilson WILSON Direct time integration for dynamic systems - Wilson-theta method [u,v,a,t] = WILSON(M,C,K,dt,p,u1,v1,[alpha delta theta]) applies the Wilson-theta method for the calculation of the nodal displacements u, velocities v and accelerations a of the dynamic system with the system matrices M, C and K due to the excitation p. M C K dt p u1 v1 u t
Mass matrix (nDof * nDof) Damping matrix (nDof * nDof) Stiffness matrix (nDof * nDof) Time step of the integration scheme (1 * 1). Should be small enough to ensure the stability and the precision of the integration scheme. Excitation (nDof * N). p(:,k) corresponds to time point t(k). Displacements at time point t(1) (nDof * 1). Defaults to zero. Velocities at time point t(1) (nDof * 1). Defaults to zero. Displacements (nDof * N). u(:,k) corresponds to time point t(k). Time axis (1 * N), defined as t = [0:N-1] * dt.
ldlt LDL Gauss factorization of symmetric matrices. [L,D] = LDL(A) returns the triangle matrix L and the vector D so that A = L.’ * diag(D) * L. You should consider renumbering A with SYMRCM before factoring the matrix. See details in K.J. Bathe, ’Finite Element Procedures in Engineering Analysis’, Prentice-Hall, 1982, p. 442
Beam functions
1.4
21
Beam functions
dof_beam ke_beam kelcs_beam trans_beam loads_beam loadslcs_beam accel_beam forces_beam forceslcs_beam nelcs_beam nedloadlcs_beam coord_beam disp_beam dispgcs2lcs_beam fdiagrgcs_beam fdiagrlcs_beam sdiagrgcs_beam sdiagrlcs_beam
Element degrees of freedom for a beam element. Beam element stiffness and mass matrix in global coordinate system. Beam element stiffness and mass matrix in local coordinate system. Transform coordinate system for a beam element. Equivalent nodal forces for a beam element in the GCS. Equivalent nodal forces for a beam element in the LCS. Compute the distributed loads for a beam due to an acceleration. Compute the element forces for a beam element. Compute the element forces for a beam element in the LCS. Shape functions for a beam element. Shape functions for a distributed load on a beam element. Coordinates of the beam elements for plotting. Return matrices to compute the displacements of the deformed beams. Transform the element displacements to the LCS for a beam. Return matrices to plot the forces in a beam element. Force diagram for a beam element in LCS. Return matrices to plot the stresses in a b eam element. Stress diagram for a beam element in LCS.
p.22 p.22 p.22 p.23 p.23 p.24 p.24 p.24 p.25 p.25 p.26 p.26 p.26 p.27 p.28 p.28 p.29 p.30
22
Functions — By category
dof beam DOF_BEAM
Element degrees of freedom for a beam element.
dof = dof_beam(NodeNum) builds the vector with the labels of the degrees of freedom for which stiffness is present in the beam element. NodeNum Node definitions dof Degrees of freedom
[NodID1 NodID2] (1 * 2) (12 * 1)
See also GETDOF.
ke beam KE_BEAM
Beam element stiffness and mass matrix in global coordinate system.
[Ke,Me] = ke_beam(Node,Section,Material,Options) returns the element stiffness and mass matrix in the global coordinate system for a two node beam element (isotropic material) Node Section Material Options Ke Me
Node definitions [x y z] (3 * 3) Section definition [A ky kz Ixx Iyy Izz] Material definition [E nu rho] Element options {Option1 Option2 ...} Element stiffness matrix (12 * 12) Element mass matrix (12 * 12)
See also KELCS_BEAM, TRANS_BEAM, ASMKM, KE_TRUSS.
kelcs beam KELCS_BEAM
Beam element stiffness and mass matrix in local coordinate system.
[KeLCS,MeLCS] = kelcs_beam(L,A,ky,kz,Ixx,Iyy,Izz,E,nu,rho,Options) [KeLCS,MeLCS] = kelcs_beam(L,A,ky,kz,Ixx,Iyy,Izz,E,nu,rho) KeLCS = kelcs_beam(L,A,ky,kz,Ixx,Iyy,Izz,E,nu) returns the element stiffness and mass matrix in the local coordinate system for a two node beam element (isotropic material) L A
Beam length Beam cross section
Beam functions
23
ky kz Ixx Iyy Izz E nu rho Options KeLCS MeLCS
Shear deflection factor A_sy = ky * A Shear deflection factor A_sz = kz * A Moment of inertia Moment of inertia Moment of inertia Young’s modulus Poisson coefficient Mass density Options for the mass matrix: {’lumped’}, {’norotaroryinertia’} Element stiffness matrix (12 * 12) Element mass matrix (12 * 12)
See also KE_BEAM, KELCS_TRUSS.
trans beam TRANS_BEAM
Transform coordinate system for a beam element.
t = trans_beam(Node) computes the transformation matrix between the local and the global coordinate system for the beam element. Node t
Node definitions Transformation matrix
[x y z] (3 * 3) (3 * 3)
See also KE_BEAM, TRANS_TRUSS.
loads beam LOADS_BEAM
Equivalent nodal forces for a beam element in the GCS.
F = loads_beam(DLoad,Node) computes the equivalent nodal forces of a distributed load (in the global coordinate system). DLoad
Distributed loads
Node F
Node definitions Load vector (12 * 1)
[n1globalX; n1globalY; n1globalZ; ...] (6 * 1) [x y z] (3 * 3)
See also LOADSLCS_BEAM, ELEMLOADS, LOADS_TRUSS.
24
Functions — By category
loadslcs beam LOADSLCS_BEAM
Equivalent nodal forces for a beam element in the LCS.
FLCS = loadslcs_beam(DLoadLCS,L) computes the equivalent nodal forces of a distributed load (in the local coordinate system). DLoadLCS
Distributed loads
L FLCS
Beam length Load vector
[n1localX; n1localY; n1localZ; ...] (6 * 1)
(12 * 1)
See also LOADS_BEAM.
accel beam ACCEL_BEAM
Compute the distributed loads for a beam due to an acceleration.
DLoads=accel_beam(Accelxyz,Elements,Sections,Materials,Options) computes the distributed loads for a beam due to an acceleration. In order to simulate gravity, accelerate the structure in the direction opposite to gravity. Accelxyz Elements Sections Materials Options DLoads
Acceleration Element definitions Section definitions Material definitions Element options Distributed loads
[Ax Ay Az] (1 * 3) [EltID TypID SecID MatID n1 n2 ...] [SecID SecProp1 SecProp2 ...] [MatID MatProp1 MatProp2 ... ] {Option1 Option2 ...} [EltID n1globalX n1globalY n1globalZ ...]
See also ACCEL, ACCEL_TRUSS.
forces beam FORCES_BEAM
Compute the element forces for a beam element.
[ForcesLCS,ForcesGCS]=forces_beam(Node,Section,Material,UeGCS,DLoad,Options) [ForcesLCS,ForcesGCS]=forces_beam(Node,Section,Material,UeGCS,DLoad) [ForcesLCS,ForcesGCS]=forces_beam(Node,Section,Material,UeGCS) computes the element forces for the beam element in the local and the global coordinate system (algebraic convention).
Beam functions
25
Node Section Material UeGCS DLoad
Node definitions Section definition Material definition Displacements (12 * 1) Distributed loads
[x y z] (3 * 3) [A ky kz Ixx Iyy Izz] [E nu]
[n1globalX; n1globalY; n1globalZ; ...] (6 * 1) Element options {Option1 Option2 ...} Element forces in the LCS (12 * 1) Element forces in the GCS (12 * 1)
Options ForcesLCS ForcesGCS
See also FORCESLCS_BEAM, ELEMFORCES.
forceslcs beam FORCESLCS_BEAM
Compute the element forces for a beam element in the LCS.
Forces=forceslcs_beam(KeLCS,UeLCS,DLoadLCS,L) Forces=forceslcs_beam(KeLCS,UeLCS) computes the element forces for the beam element in the local coordinate system (algebraic convention). KeLCS UeLCS DLoadLCS Forces
Element stiffness matrix (12 * 12) Displacements (12 * 1) Distributed loads [n1localX; n1localY; n1localZ; ...] Element forces [N; Vy; Vz; T; My; Mz](12 * 1)
See also FORCES_BEAM, FORCESLCS_TRUSS
nelcs beam NELCS_BEAM
Shape functions for a beam element.
NeLCS = nelcs_beam(Points) NeLCS = nelcs_beam(Points,phi_y,phi_z) determines the values of the shape functions in the specified points. Points NeLCS phi_y phi_z
Points in the local coordinate system (1 * nPoints) Values (nPoints * 12) Shear deformation constant in y dir (scalar) Shear deformation constant in z dir (scalar)
26
Functions — By category
See also DISP_BEAM.
nedloadlcs beam NEDLOADLCS_BEAM Shape functions for a distributed load on a beam element. NeLCS = nedloadlcs_beam(Points) NeLCS = nedloadlcs_beam(Points,phi_y,phi_z) determines the values of the shape functions for a distributed load in the specified points. These are used to compute the displacements that occur due to the distributed loads if all nodes are fixed. Points NeDLoad phi_y phi_z
Points in the local coordinate system (1 * nPoints) Values (nPoints * 6) Shear deformation constant in y dir (scalar) Shear deformation constant in z dir (scalar)
See also DISP_BEAM, NELCS_BEAM.
coord beam COORD_BEAM
Coordinates of the beam elements for plotting.
coord_beam(Nodes,NodeNum) returns the coordinates of the beam elements for plotting. Nodes Node definitions [NodID x y z] (nNodes * 4) NodeNum Node numbers [NodID1 NodID2 NodID3] (nElem * 3) X X coordinates (2 * nElem) Y Y coordinates (2 * nElem) Z Z coordinates (2 * nElem) See also COORD_TRUSS, PLOTELEM.
disp beam DISP_BEAM
Return matrices to compute the displacements of the deformed beams.
Beam functions
27
[Ax,Ay,Az,B,Cx,Cy,Cz] =disp_beam(Nodes,Elements,DOF,EltIDDLoad,Sections,Materials,Points) [Ax,Ay,Az,B,Cx,Cy,Cz] =disp_beam(Nodes,Elements,DOF,EltIDDLoad,Sections,Materials) [Ax,Ay,Az,B] =disp_beam(Nodes,Elements,DOF,[],Sections,Materials) [Ax,Ay,Az,B] =disp_beam(Nodes,Elements,DOF) returns the matrices to compute the displacements of the deformed beams. The coordinates of the specified points along the deformed beams element are computed using X=Ax*U+Cx*DLoad+B(:,1); Y=Ay*U+Cy*DLoad+B(:,2) and Z=Az*U+Cz*DLoad+B(:,3). The matrices Cx,Cy and Cz superimpose the displacements that occur due to the distributed loads if all nodes are fixed. Nodes Elements DOF EltIDDLoad Sections Materials Points Ax Ay Az B Cx Cy Cz
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Degrees of freedom (nDOF * 1) Elements with distributed loads [EltID] (use empty array [] when shear deformation is considered but no DLoads are pres Section definitions [SecID SecProp1 SecProp2 ...] Material definitions [MatID MatProp1 MatProp2 ... ] Points in the local coordinate system (1 * nPoints) Matrix to compute the x-coordinates of the deformations Matrix to compute the y-coordinates of the deformations Matrix to compute the z-coordinates of the deformations Matrix which contains the x-, y- and z-coordinates of the undeformed structure Matrix to compute the x-coordinates of the deformations Matrix to compute the y-coordinates of the deformations Matrix to compute the z-coordinates of the deformations
See also DISP_TRUSS, PLOTDISP, NELCS_BEAM, NEDLOADLCS_BEAM.
dispgcs2lcs beam DISPGCS2LCS_BEAM
Transform the element displacements to the LCS for a beam.
UeLCS=dispgcs2lcs_beam(UeGCS,Node) transforms the element displacements from the GCS to the LCS for a beam element. Node UeGCS UeLCS
Node definitions [x y z] (3 * 3) Displacements in the GCS (12 * 1) Displacements in the LCS (12 * 1)
28
Functions — By category
See also DISPGCS2LCS_TRUSS.
fdiagrgcs beam FDIAGRGCS_BEAM
Return matrices to plot the forces in a beam element.
[ElemGCS,FdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = fdiagrgcs_beam(ftype,Forces,Node,DLoad,Points) [ElemGCS,FdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = fdiagrgcs_beam(ftype,Forces,Node,DLoad) returns the coordinates of the points along the beam in the global coordinate system and the coordinates of the forces with respect to the beam in the global coordinate system. These can be added in order to plot the forces: ElemGCS+FdiagrGCS. The coordinates of the points with extreme values and the coordinates of the extreme values with respect to the beam are given as well and can be similarly added: ElemExtGCS+ExtremaGCS. Extrema is the list with the correspondig extreme values. ftype
’norm’ Normal force (in the local x-direction) ’sheary’ Shear force in the local y-direction ’shearz’ Shear force in the local z-direction ’momx’ Torsional moment (around the local x-direction) ’momy’ Bending moment around the local y-direction ’momz’ Bending moment around the local z-direction Forces Element forces in LCS (beam convention) [N; Vy; Vz; T; My; Mz] (12 * 1) Node Node definitions [x y z] (3 * 3) DLoad Distributed loads [n1globalX; n1globalY; n1globalZ; ...] (6 * 1) Points Points in the local coordinate system (1 * nPoints) ElemGCS Coordinates of the points along the beam in GCS (nPoints * 3) FdiagrGCS Coordinates of the force with respect to the beam in GCS (nValues * 3) ElemExtGCS Coordinates of the points with extreme values in GCS (nValues * 3) ExtremaGCS Coordinates of the extreme values with respect to the beam in GCS (nValues * 3) Extrema Extreme values (nValues * 1) See also PLOTFORC, FDIAGRLCS_BEAM, FDIAGRGCS_TRUSS.
fdiagrlcs beam FDIAGRLCS_BEAM
Force diagram for a beam element in LCS.
Beam functions
29
[FdiagrLCS,loc,Extrema] = fdiagrlcs_beam(ftype,Forces,DLoadLCS,L,Points) computes the elements forces at the specified points. The extreme values are analytically determined. ftype
Forces DLoadLCS Points FdiagrLCS loc Extrema
’norm’ Normal force (in the local x-direction) ’sheary’ Shear force in the local y-direction ’shearz’ Shear force in the local z-direction ’momx’ Torsional moment (around the local x-direction) ’momy’ Bending moment around the local y-direction ’momz’ Bending moment around the local z-direction Element forces in LCS (beam convention) [N; Vy; Vz; T; My; Mz](12 * 1) Distributed loads in LCS [n1localX; n1localY; n1localZ; ...](6 * 1) Points in the local coordinate system (1 * nPoints) Element forces at the points (1 * nPoints) Locations of the extreme values (nValues * 1) Extreme values (nValues * 1)
See also FDIAGRGCS_BEAM.
sdiagrgcs beam SDIAGRGCS_BEAM
Return matrices to plot the stresses in a beam element.
[ElemGCS,SdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = sdiagrgcs_beam(stype,Forces,Node,DLoad,Section,Points) [ElemGCS,SdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = sdiagrgcs_beam(stype,Forces,Node,DLoad,Section) returns the coordinates of the points along the beam in the global coordinate system and the coordinates of the stresses with respect to the beam in the global coordinate system. These can be added in order to plot the stresses: ElemGCS+SdiagrGCS. The coordinates of the points with extreme values and the coordinates of the extreme values with respect to the beam are given as well and can be similarly added: ElemExtGCS+ExtremaGCS. Extrema is the list with the correspondig extreme values. stype
’snorm’ ’smomyt’ ’smomyb’ ’smomzt’ ’smomzb’ ’smax’
Normal stress due to normal force Normal stress due to bending moment around the local y-direction at the Normal stress due to bending moment around the local y-direction at the Normal stress due to bending moment around the local z-direction at the Normal stress due to bending moment around the local z-direction at the Maximal normal stress (normal force
top bottom top bottom and bending moment)
30
Functions — By category
’smin’ Minimal normal stress (normal force and bending moment) Forces Element forces in LCS (beam convention) [N; Vy; Vz; T; My; Mz] (12 * 1) Node Node definitions [x y z] (3 * 3) DLoad Distributed loads [n1globalX; n1globalY; n1globalZ; ...] (6 * 1) Section Section definition [A ky kz Ixx Iyy Izz yt yb zt zb] Points Points in the local coordinate system (1 * nPoints) ElemGCS Coordinates of the points along the beam in GCS (nPoints * 3) SdiagrGCS Coordinates of the stress with respect to the beam in GCS (nValues * 3) ElemExtGCS Coordinates of the points with extreme values in GCS (nValues * 3) ExtremaGCS Coordinates of the extreme values with respect to the beam in GCS (nValues * 3) Extrema Extreme values (nValues * 1) See also PLOTSTRESS, SDIAGRLCS_BEAM, SDIAGRGCS_TRUSS.
sdiagrlcs beam SDIAGRLCS_BEAM
Stress diagram for a beam element in LCS.
[SdiagrLCS,loc,Extrema] = sdiagrlcs_beam(ftype,Forces,DLoadLCS,L,Points) computes the stresses at the specified points. The extreme values are analytically determined. stype
Forces DLoadLCS Points Section SdiagrLCS loc Extrema
’snorm’ ’smomyt’
Normal stress due to normal force Normal stress due to bending moment around the local y-direction at the top ’smomyb’ Normal stress due to bending moment around the local y-direction at the bottom ’smomzt’ Normal stress due to bending moment around the local z-direction at the top ’smomzb’ Normal stress due to bending moment around the local z-direction at the bottom ’smax’ Maximal normal stress (normal force and bending moment) ’smin’ Minimal normal stress (normal force and bending moment) Element forces in LCS (beam convention) [N; Vy; Vz; T; My; Mz](12 * 1) Distributed loads in LCS [n1localX; n1localY; n1localZ; ...](6 * 1) Points in the local coordinate system (1 * nPoints) Section definition [A ky kz Ixx Iyy Izz yt yb zt zb] Stresses at the points (1 * nPoints) Locations of the extreme values (nValues * 1) Extreme values (nValues * 1)
See also SDIAGRGCS_BEAM.
Truss functions
1.5
31
Truss functions
dof_truss Element degrees of freedom for a truss element. Truss element stiffness and mass matrix in global coordinate system. ke_truss kelcs_truss Truss element stiffness and mass matrix in local coordinate system. trans_truss Transform coordinate system for a truss element. loads_truss Equivalent nodal forces for a truss element in the GCS. accel_truss Compute the distributed loads for a truss due to an acceleration. Compute the element forces for a truss element. forces_truss forceslcs_truss Compute the element forces for a truss element in the LCS. coord_truss Coordinates of the truss elements for plotting. disp_truss Return matrices to compute the displacements of the deformed trusses. dispgcs2lcs_trussTransform the element displacements to the LCS for a truss. fdiagrgcs_truss Return matrices to plot the forces in a truss element. sdiagrgcs_truss Return matrices to plot the stresses in a truss element.
p.32 p.32 p.32 p.33 p.33 p.33 p.34 p.34 p.35 p.35 p.36 p.36 p.37
32
Functions — By category
dof truss DOF_TRUSS
Element degrees of freedom for a truss element.
dof = dof_truss(NodeNum) builds the vector with the labels of the degrees of freedom for which stiffness is present in the truss element. NodeNum Node numbers [NodID1 NodID2] (1 * 2) dof Degrees of freedom (6 * 1) See also GETDOF.
ke truss KE_TRUSS
Truss element stiffness and mass matrix in global coordinate system.
[Ke,Me] = ke_truss(Node,Section,Material,Options) returns the element stiffness and mass matrix in the global coordinate system for a two node truss element (isotropic material) Node Section Material Options Ke Me
Node definitions Section definition Material definition Element options Element stiffness matrix Element mass matrix
[x1 y1 z1; x2 y2 z2] (2 * 3) [A] [E nu rho] {Option1 Option2 ...} (6 * 6) (6 * 6)
See also KELCS_TRUSS, TRANS_TRUSS, ASMKM, KE_BEAM.
kelcs truss KELCS_TRUSS
Truss element stiffness and mass matrix in local coordinate system.
[KeLCS,MeLCS] = kelcs_truss(L,A,E,rho,Options) KeLCS = kelcs_truss(L,A,E) returns the element stiffness and mass matrix in the local coordinate system for a two node truss element (isotropic material) L A E
Truss length Truss cross section Young’s modulus
Truss functions
rho Options KeLCS MeLCS
33
Mass density Options for the mass matrix: {’lumped’} Element stiffness matrix (6 * 6) Element mass matrix (6 * 6)
See also KE_TRUSS, KELCS_BEAM.
trans truss TRANS_TRUSS
Transform coordinate system for a truss element.
t = trans_truss(Node) computes the transformation matrix between the local and the global coordinate system for the truss element. Node t
Node definitions Transformation matrix
[x y z] (2 * 3) (3 * 3)
See also KE_TRUSS, TRANS_BEAM.
loads truss LOADS_TRUSS
Equivalent nodal forces for a truss element in the GCS.
F = loads_truss(DLoad,Node) computes the equivalent nodal forces of a distributed load (in the global coordinate system). DLoad
Distributed loads
Node F
Node definitions Load vector (6 * 1)
[n1globalX n1globalY n1globalZ ...] (6 * 1) [x y z] (2 * 3)
See also ELEMLOADS, LOADS_BEAM.
accel truss ACCEL_TRUSS
Compute the distributed loads for a truss due to an acceleration.
34
Functions — By category
DLoads=accel_truss(Accelxyz,Elements,Sections,Materials,Options) computes the distributed loads for a truss due to an acceleration. In order to simulate gravity, accelerate the structure in the direction opposite to gravity. Accelxyz Elements Sections Materials Options DLoads
Acceleration Element definitions Section definitions Material definitions Element options Distributed loads
[Ax Ay Az] (1 * 3) [EltID TypID SecID MatID n1 n2 ...] [SecID SecProp1 SecProp2 ...] [MatID MatProp1 MatProp2 ... ] {Option1 Option2 ...} [EltID n1globalX n1globalY n1globalZ ...]
See also ACCEL, ACCEL_BEAM.
forces truss FORCES_TRUSS
Compute the element forces for a truss element.
[ForcesLCS,ForcesGCS]=forces_truss(Node,Section,Material,UeGCS) computes the element forces for the truss element in the local and the global coordinate system (algebraic convention). Node Section Material UeGCS Options ForcesLCS ForcesGCS
Node definitions [x y z] (2 * 3) Section definition [A] Material definition [E] Displacements (6 * 1) Element options {Option1 Option2 ...} Element forces in the LCS (12 * 1) Element forces in the GCS (12 * 1)
See also FORCESLCS_TRUSS, ELEMFORCES.
forceslcs truss FORCESLCS_TRUSS
Compute the element forces for a truss element in the LCS.
Forces=forceslcs_truss(KeLCS,UeLCS) computes the element forces for the truss element in the local coordinate system (algebraic convention). KeLCS
Element stiffness matrix (6 * 6)
Truss functions
UeLCS Forces
35
Displacements (6 * 1) Element forces
[N; 0; 0](6 * 1)
See also FORCES_TRUSS, FORCESLCS_BEAM
coord truss COORD_TRUSS
Coordinates of the truss elements for plotting.
coord_truss(Nodes,NodeNum) returns the coordinates of the truss elements for plotting. Nodes Node definitions [NodID x y z] (nNodes * 4) NodeNum Node numbers [NodID1 NodID2] (nElem * 2) X X coordinates (2 * nElem) Y Y coordinates (2 * nElem) Z Z coordinates (2 * nElem) See also COORD_BEAM, PLOTELEM.
disp truss DISP_TRUSS
Return matrices to compute the displacements of the deformed trusses.
[Ax,Ay,Az,B]=disp_truss(Nodes,Elements,DOF,[],[],[],Points) [Ax,Ay,Az,B]=disp_truss(Nodes,Elements,DOF) returns the matrices to compute the displacements of the deformed trusses. The coordinates of the specified points along the deformed beams element are computed using X=Ax*U+B(:,1); Y=Ay*U+B(:,2) and Z=Az*U+B(:,3). Nodes Elements DOF Points Ax Ay Az B
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Degrees of freedom (nDOF * 1) Points in the local coordinate system (1 * nPoints) Matrix to compute the x-coordinates of the deformations Matrix to compute the y-coordinates of the deformations Matrix to compute the z-coordinates of the deformations Matrix which contains the x-, y- and z-coordinates of the undeformed structure
See also DISP_BEAM, PLOTDISP.
36
Functions — By category
dispgcs2lcs truss DISPGCS2LCS_TRUSS
Transform the element displacements to the LCS for a truss.
UeLCS=dispgcs2lcs_truss(UeGCS,Node) transforms the element displacements from the GCS to the LCS for a truss element. Node UeGCS UeLCS
Node definitions [x y z] (2 * 3) Displacements in the GCS (6 * 1) Displacements in the LCS (6 * 1)
See also DISPGCS2LCS_BEAM.
fdiagrgcs truss FDIAGRGCS_TRUSS
Return matrices to plot the forces in a truss element.
[ElemGCS,FdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = fdiagrgcs_truss(ftype,Forces,Node,[],Points) [ElemGCS,FdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = fdiagrgcs_truss(ftype,Forces,Node) returns the coordinates of the points along the truss in the global coordinate system and the coordinates of the forces with respect to the truss in the global coordinate system. These can be added in order to plot the forces: ElemGCS+FdiagrGCS. The coordinates of the points with extreme values and the coordinates of the extreme values with respect to the truss are given as well and can be similarly added: ElemExtGCS+ExtremaGCS. Extrema is the list with the correspondig extreme values. ftype Forces Node Points ElemGCS FdiagrGCS
’norm’ Normal force (in the local x-direction) Element forces in LCS [N; 0; 0; 0; 0; 0](12 * 1) Node definitions [x y z] (3 * 3) Points in the local coordinate system (1 * nPoints) Coordinates of the points along the truss in GCS (nPoints * 3) Coordinates of the force with respect to the truss in GCS (nValues * 3) ElemExtGCS Coordinates of the points with extreme values in GCS (nValues * 3) ExtremaGCS Coordinates of the extreme values with respect to the truss in GCS (nValues * 3) Extrema Extreme values (nValues * 1) See also PLOTFORC, FDIAGRGCS_BEAM.
Truss functions
37
sdiagrgcs truss SDIAGRGCS_TRUSS
Return matrices to plot the stresses in a truss element.
[ElemGCS,SdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = sdiagrgcs_truss(stype,Forces,Node,[],Section,Points) [ElemGCS,SdiagrGCS,ElemExtGCS,ExtremaGCS,Extrema] = sdiagrgcs_truss(stype,Forces,Node,[],Section) returns the coordinates of the points along the truss in the global coordinate system and the coordinates of the stresses with respect to the truss in the global coordinate system. These can be added in order to plot the stresses: ElemGCS+SdiagrGCS. The coordinates of the points with extreme values and the coordinates of the extreme values with respect to the truss are given as well and can be similarly added: ElemExtGCS+ExtremaGCS. Extrema is the list with the correspondig extreme values. stype Forces Node Section Points ElemGCS SdiagrGCS
’snorm’ Normal stress due to normal force Element forces in LCS [N; 0; 0; 0; 0; 0](12 * 1) Node definitions [x y z] (3 * 3) Section definition [A ky kz Ixx Iyy Izz yt yb zt zb] Points in the local coordinate system (1 * nPoints) Coordinates of the points along the truss in GCS (nPoints * 3) Coordinates of the force with respect to the truss in GCS (nValues * 3) ElemExtGCS Coordinates of the points with extreme values in GCS (nValues * 3) ExtremaGCS Coordinates of the extreme values with respect to the truss in GCS (nValues * 3) Extrema Extreme values (nValues * 1) See also PLOTSTRESS, SDIAGRGCS_BEAM.
38
1.6
Functions — By category
General shell functions
elempressure Equivalent nodal forces for pressure ⊥ shell surface. Compute the shell forces and moments (element solution). elemshellf elemstress Compute the element stresses (element solution). gaussq Gauss points for 2D numerical integration. makemesh Creates a mesh of quadrilateral elements for a surface defined by 4 lines. nodalshellf Compute the nodal shell forces/moments from element solution. Compute the nodal stresses from element solution. nodalstress plotprincstress Plot the principal stresses. plotshellfcontour Plot contour lines of forces/moments per unit length. plotstresscontour Plot stress contour lines. plotshellfcontourfPlot filled contours of shell forces. Compute the principal stresses and directions. principalstress printshellf Display forces/moments in command window (shell elements). printstress Display stresses in command window (shell elements).
p.39 p.39 p.39 p.40 p.41 p.42 p.42 p.43 p.43 p.44 p.44 p.45 p.46 p.46
General shell functions
39
elempressure ELEMPRESSURE
Equivalent nodal forces for a pressure load on a shell element.
F = elempressure(Pressures,Nodes,Elements,Types,DOF) computes the equivalent nodal forces of a distributed pressure (in the local coordinate system xyz, with z perpendicular to the surface). Pressures Pressure perpendicular to surface or edge Nodes Node definitions Elements Element definitions Types Element type definitions DOF Degrees of freedom F Load vector
[EltID n1localZ n2localZ n3localZ .. [NodID x y z] [EltID TypID SecID MatID n1 n2 ...] {TypID EltName Option1 ... } (nDOF * 1) (nDOF * 1)
See also PRESSURE_SHELL8, PRESSURE_SHELL4, NODALVALUES.
elemshellf ELEMSHELLF
Compute the shell forces and moments per unit length for shell elements (element solution).
[FeLCS] = elemshellf(Elements,Sections,SeLCS) computes the element forces/moments per unit length in the local coordinate system. Elements Sections SeLCS
Element definitions [EltID TypID SecID MatID n1 n2 ...] Section definitions [SecID SecProp1 SecProp2 ...] Element stresses in LCS in corner nodes IJKL and at top/mid/bot of shell (nElem * 72) [sxx syy szz sxy syz sxz] FeLCS Element forces/moments per unit length (nElem * 32) (32 = 8 forces * 4 corner nodes) [Nx Ny Nxy Mx My Mxy Vx Vy] See also ELEMSTRESS.
elemstress ELEMSTRESS
Compute the element stresses.
[SeGCS,SeLCS,vLCS] = elemstress(Nodes,Elements,Types,Sections,Materials,DOF,U) [SeGCS,SeLCS] = elemstress(Nodes,Elements,Types,Sections,Materials,DOF,U) SeGCS = elemstress(Nodes,Elements,Types,Sections,Materials,DOF,U)
40
Functions — By category
computes the element stresses in the global and the local coordinate system. Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } Sections Section definitions [SecID SecProp1 SecProp2 ...] Materials Material definitions [MatID MatProp1 MatProp2 ... ] DOF Degrees of freedom (nDOF * 1) U Displacements (nDOF * 1) elemstress(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’GCS’ Sets the gcs in which SeGCS is returned. Default: ’cart’. Type of values: ’cart’: cartesian coordinate system ’cyl’ : cylindrical coordinate system ’sph’ : spherical coordinate system SeGCS
SeLCS
vLCS
Element stresses in GCS in corner nodes IJKL and at top/mid/bot of shell (nElem * 72) 72 = 6 stresscomp. * 4 nodes * 3 locations [sxx syy szz sxy syz sxz] Element stresses in LCS in corner nodes IJKL and at top/mid/bot of shell (nElem * 72) [sxx syy szz sxy syz sxz] Unit vectors of LCS (nElem * 9)
See also SE_SHELL8, SE_SHELL4.
gaussq GAUSSQ
Gauss points for 2D numerical integration.
[x,H] = gaussq(n) returns the coordinates and weights for a 2D gauss-legendre quadrature. n x H
number of points in one direction = 2 or 3 coordinates of gauss-points (n^2 * 2) weights used in summation (1 * n^2)
See also KE_SHELL8, KELCS_SHELL4
General shell functions
41
makemesh MAKEMESH
Creates a mesh of quadrilateral elements for a surface defined by 4 lines
[Nodes,Elements,Edge1,Edge2,Edge3,Edge4] = makemesh(Line1,Line2,Line3,Line4,m,n,Type,Section,Material) Line1,2,3,4 Lines that define the surface. (n1 * 3) These lines: -should be defined in a clockwise or counter-clockwise order (this means: first point of Line2 is the last point of Line1 in the case of a counter-clockwise direction) -can be defined by any number of points m n Type Section Material
Number of elements that divide Line2 and Line4 Number of elements that divide Line1 and Line3 Cell containing Type number and name {TypeID EltName} Section ID Material ID
Nodes Elements Edge1,2,3,4 Normals
Nodes definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] vector containing Nodes on Line1,2,.. ((m+1) * 1) or ((n+1) * 1) Normals of shell surface in Nodes [NodID nx ny nz]
makemesh(...,ParamName,ParamValue)sets the value of the specified parameters. The following parameters can be specified: ’L1method’ Specify the interpolation scheme for Line1. default: ’spline’. This method is used in the MATLAB INTERP1 function. See DOC INTERP1 for available interpolation schemes. ’L2method’ ’L3method’ ’L4method’ ’UniformSpacing’ {’on’,default: ’off’} If UniformSpacing is set to ’off’, the distance between consecutive nodes along the edges of the mesh is based on the spacing between the points in the corresponding Line. If UniformSpacing is ’on’, makemesh attempts to distribute the nodes uniformly along the lines regardless of the initial spacing between the points by approximating the natural parametrizations of the lines.
See also GRID_SHELL8, GRID_SHELL4.
42
Functions — By category
nodalshellf NODALSHELLF
Compute the nodal shell forces/moments per unit length from the element solutio
[FnLCS,FnLCS2] = nodalshellf(Nodes,Elements,Types,FeLCS) FnLCS = nodalshellf(Nodes,Elements,Types,FeLCS) computes the nodal forces from the element solution. Nodes Elements Types Sections FeLCS
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Section definitions [SecID SecProp1 SecProp2 ...] Element forces/moments per unit length in corner nodes IJKL (nElem * 32) [Nx Ny Nxy Mx My Mxy Vx Vy] FnLCS Nodal forces/moments per unit length in corner nodes IJKL (nElem * 32) [Nx Ny Nxy Mx My Mxy Vx Vy] FnLCS2 Nodal forces per node (nNodes * 9) (9 = NodID + 8 fcomp.) [NodID Nx Ny Nxy Mx My Mxy Vx Vy] See also ELEMSHELLF.
nodalstress NODALSTRESS
Compute the nodal stresses from the element solution.
[Sn,Sn2] = nodalstress(Nodes,Elements,Types,Se) Sn = nodalstress(Nodes,Elements,Types,Se) computes the nodal stresses from the element solution. Nodes Elements Types Sections Se
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Section definitions [SecID SecProp1 SecProp2 ...] Element stresses in corner nodes IJKL and at top/mid/bot of shell (nElem * 72) 72 = 6 stresscomp. * 4 nodes * 3 locations [sxx syy szz sxy syz sxz ...] Sn Nodal stress in corner nodes IJKL and at top/mid/bot of shell (nElem * 72) 72 = 6 stresscomp. * 4 nodes * 3 locations [sxx syy szz sxy syz sxz ...] Sn2 Nodal stresses per node at top/mid/bot of shell (nNodes * 19) (19 = NodID + 3 locations * 6 scomp.) [NodID sxx syy szz sxy syz sxz ...] See also ELEMSTRESS.
General shell functions
43
plotprincstress PLOTPRINCSTRESS
Plot the principal stresses in shell elements.
plotprincstress(Nodes,Elements,Types,Spr,Vpr) plots the principal stresses in shell elements with a vector plot. Nodes Elements Types Spr Vpr
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Element type definitions {TypID EltName Option1 ... } Principal stresses (nElem * 72) [s3 s2 s1 0 0 0 ...] Principal dir. matrix (output from principalstress) {nElem * 12} plotprincdir(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’location’ Location (top,mid,bot). Default: ’top’. ’GCS’ Plot the GCS. Default: ’on’. ’VectScal’ Vector scaling. Default: ’auto’. ’Undeformed’ Plots the undeformed mesh. Default: ’k-’. ’Handle’ Plots in the axis with this handle. Default: current axis. See also PRINCIPALSTRESS.
plotshellfcontour PLOTSHELLFCONTOUR
Plot forces/moments per unit length in shell elements.
plotshellfcontour(ftype,Nodes,Elements,Types,F) plots force contours (output from ELEMSHELLF/NODALSHELLF). ftype
’Nx’ Membrane forces ’Ny’ ’Nxy’ ’Mx’ Bending moments ’My’ ’Mxy’ ’Vx’ Shear forces ’Vy’ Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } F Forces/moments per unit length (nElem * 32) [Nx Ny Nxy Mx My Mxy Vx Vy] plotshellfcontour(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified:
44
Functions — By category
’Ncontour’ ’GCS’ ’Undeformed’ ’Handle’
Number of contours. Default: ’10’. Plot the GCS. Default: ’on’. Plots the undeformed mesh. Default: ’k-’. Plots in the axis with this handle. Default: current axis.
See also ELEMSHELLF, SCONTOUR_SHELL8, SCONTOUR_SHELL4.
plotstresscontour PLOTSTRESSCONTOUR
Plot stress contour lines in shell elements.
plotstresscontour(stype,Nodes,Elements,Types,S) plots stress contours (output from ELEMSTRESS/NODALSTRESS). stype
’sx’ Normal stress (in the global/local x-direction) ’sy’ ’sz’ ’sxy’ Shear stress ’syz’ ’sxz’ Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } S Element stresses in GCS/LCS [sxx syy szz sxy syz sxz] (nElem * 72) plotstresscontour(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’location’ Location (top,mid,bot). Default: ’top’. ’Ncontour’ Number of contours. Default: ’10’. ’GCS’ Plot the GCS. Default: ’on’. ’Undeformed’ Plots the undeformed mesh. Default: ’k-’. ’Handle’ Plots in the axis with this handle. Default: current axis. See also ELEMSTRESS, SCONTOUR_SHELL8, SCONTOUR_SHELL4.
plotshellfcontourf PLOTSHELLFCONTOURF
Plot filled contours of shell forces in shell elements.
plotshellfcontourf(stype,Nodes,Elements,Types,F) plotshellfcontourf(stype,Nodes,Elements,Types,F,DOF,U) plots force contours (output from ELEMSTRESS).
General shell functions
45
ftype
’Nx’ Membrane forces ’Ny’ ’Nxy’ ’Mx’ Bending moments ’My’ ’Mxy’ ’Vx’ Shear forces ’Vy’ Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } F Forces/moments per unit length (nElem * 32) [Nx Ny Nxy Mx My Mxy Vx Vy] DOF Degrees of freedom (nDOF * 1) U Displacements (nDOF * 1) plotshellfcontourf(...,ParamName,ParamValue) sets the value of the specified parameters. The following parameters can be specified: ’GCS’ Plot the GCS. Default: ’on’. ’ncolor’ Number of colors in colormap. Default: 10. ’Handle’ Plots in the axis with this handle. Default: current axis. See also ELEMSHELLF, PATCH_SHELL8, PATCH_SHELL4.
principalstress PRINCIPALSTRESS
Compute the principal stresses and directions in shell elements.
[Spr,Vpr] = principalstress(Elements,SeGCS) computes the principal stresses and directions. Elements SeGCS
Spr Vpr
Element definitions [EltID TypID SecID MatID n1 n2...] Element stresses in GCS in corner nodes IJKL and at top/mid/bot of shell (nElem * 72) 72 = 6 stresscomp. * 4 nodes * 3 locations [sxx syy szz sxy syz sxz ...] Principal stress (nElem * 72) [s1 s2 s3 0 0 0 ...] Principal stress directions {nElem * 12}
See also ELEMSTRESS, PLOTPRINCDIR.
46
Functions — By category
printshellf PRINTSHELLF Display forces/moments in command window (shell elements). printshellf(Elements,F) displays shell forces in command window. Elements F
Element definitions Shellf matrix
[EltID TypID SecID MatID n1 n2 ...] [Nx Ny Nxy Mx My Mxy Vx Vy] (nElem * 32)
See also PRINTSTRESS.
printstress PRINTSTRESS Display stress in command window (shell elements). printstress(Elements,S,location) printstress(Elements,S) displays stress in command window. Elements S
Element definitions Stress matrix
See also PRINTSHELLF.
[EltID TypID SecID MatID n1 n2 ...] [sx sy sz sxy syz szx] (nElem * 72)
Shell4 functions
1.7
47
Shell4 functions
dof_shell4 Element degrees of freedom for a shell4 element. Shell4 element stiffness and mass matrix in global coordinate system. ke_shell4 kelcs_shell4 Shell4 element stiffness and mass matrix in local coordinate system. trans_shell4 Transform coordinate system for a shell4 element. ke_dkt DKT plate element stiffness and mass matrix. q_dkt Q matrix for a DKT element. Shape functions for a quadrilateral serendipity element with 4 nodes. sh_qs4 sh_t Shape functions for a triangular plate element. se_shell4 Compute the element stresses for a shell4 element. selcs_shell4 Compute the element stresses for a shell4 element. accel_shell4 Compute the distributed loads for a shell due to an acceleration. Equivalent nodal forces for a shell4 element in the GCS. loads_shell4 loadslcs_shell4 Equivalent nodal forces for a shell4 element in the LCS. pressure_shell4 Equivalent nodal forces for a shell4 element in the GCS due to a pressure. coord_shell4 Coordinates of the shell4 elements for plotting. Matrices to compute the displacements of the deformed shell. disp_shell4 scontour_shell4 Matrix to plot contours in a shell4 element. patch_shell4 Patch information of the shell4 elements for plotting. grid_shell4 Grid in natural coordinates for mapped meshing.
p.48 p.48 p.49 p.49 p.50 p.50 p.51 p.51 p.53 p.54 p.51 p.52 p.52 p.53 p.54 p.55 p.55 p.56 p.56
48
Functions — By category
dof shell4 DOF_SHELL4
Element degrees of freedom for a shell4 element.
dof = dof_shell4(NodeNum) builds the vector with the labels of the degrees of freedom for which stiffness is present in the shell4 element. NodeNum Node definitions dof Degrees of freedom
[NodID1 NodID2 ... NodIDn]
(1 * 4) (24 * 1)
See also GETDOF.
ke shell4 KE_SHELL4
shell element stiffness and mass matrix in global coordinate system.
[Ke,Me] = ke_shell4(Node,Section,Material,Options) Ke = ke_shell4(Node,Section,Material,Options) returns the element stiffness and mass matrix in the global coordinate system for a four node shell element (isotropic material). Node
Node definitions [x y z] (4 * 3) Nodes should have the following order: 4---------3 | | | | | | 1---------2
Section Material Options Ke Me
Section definition [h] Material definition [E nu rho] Element options {Option1 Option2 ...} Element stiffness matrix (24 * 24) Element mass matrix (24 * 24)
This shell element consists of a bilinear membrane element and four overlaid DKT triangles for the bending stiffness. See also KE_BEAM, ASMKM, KE_TRUSS.
Shell4 functions
49
kelcs shell4 KELCS_SHELL4
shell element stiffness and mass matrix in element coordinate system.
[Ke,Me] = kelcs_shell4(Node_lc,h,E,nu,rho) Ke = kelcs_shell4(Node_lc,h,E,nu) returns the element stiffness and mass matrix in the element coordinate system for a four node shell element (isotropic material). Node
Node definitions [x y z] (4 * 3) Nodes should have the following order: 4---------3 | | | | | | 1---------2
h E nu rho Options KeLCS MeLCS
Shell thickness Young’s modulus Poisson coefficient Mass density Element options [NOT available yet] Element stiffness matrix (24 * 24) Element mass matrix (24 * 24)
{Option1 Option2 ...}
This element is a flat shell element that consists of a bilinear membrane element and four overlaid DKT triangles for the bending stiffness. See also KE_SHELL8, ASMKM, KE_DKT.
trans shell4 TRANS_SHELL4
Transform coordinate system for a shell4 element.
[t,Node_lc,W] = trans_shell4(Node) [t,Node_lc] = trans_shell4(Node) t = trans_shell4(Node) computes the transformation matrix between the local and the global coordinate system and the correction matrix for non-coplanar nodes for the shell4 element. Node t Node_lc
Node definitions Transformation matrix Nodes in LCS
[x y z] (4 * 3) (3 * 3) [x y z] (4 * 3)
50
Functions — By category
W
Correction matrix for warped elements (24 * 24)
See also KE_BEAM, TRANS_TRUSS.
ke dkt KE_DKT
DKT plate element stiffness and mass matrix.
[Ke,Me] = ke_dkt(Node,h,E,nu,rho) Ke = ke_dkt(Node,h,E,nu) returns the element stiffness and mass matrix in the global coordinate system for a three node plate element (isotropic material) in the xy-plane. Node h E nu rho Ke Me
Node definitions [x y] (3 * 2) Plate thickness (uniform or defined in nodes) [h]/[h1 h2 h3] Young’s modulus Poisson coefficient Mass density Element stiffness matrix (24 * 24) Element mass matrix (24 * 24)
This element is a Discrete Kirchoff Triangle plate element. This code is a MATLAB code based on the E-2 element FORTRAN coding, described in: Construction of new efficient three-node triangular thin plate bending elements, C. Jeyachandrabose and J. Kirkhope, Computers & Structures Vol. 23, No. 5, pp. 587-603, 1986. See also Q_DKT, SH_T, KE_DKT4.
q dkt Q_DKT
Q matrix for a DKT element.
Q = Q_DKT(b,c,det) returns the Q matrix of the DKT element, which is used to compute the curvatures of the plate element b c det
Geometrical property of the triangle (see ke_dkt) (3 * 1) Geometrical property of the triangle (3 * 1) Determinant of the parametric transformation (=2*area of triangle)
Shell4 functions
51
See also KE_DKT, KELCS_SHELL4, SE_SHELL4.
sh qs4 SH_QS4
Shape functions for a quadrilateral serendipity element with 4 nodes.
[Ni,dN_dxi,dN_deta] = sh_qs4(xi,eta) returns the shape functions and its derivatives to the natural coordinates in the point (xi,eta). xi eta Ni dN_dxi dN_deta
Natural coordinate Natural coordinate Shape functions in point (xi,eta) Derivative of Ni to xi Derivative of Ni to eta
(scalar) (scalar) (4 * 1) (4 * 1) (4 * 1)
See also KELCS_SHELL4.
sh t SH_T
Shape functions for a triangular plate element.
[Ni] = sh_t(L,b,c) returns the shape functions and its derivatives in point L. L b c Ni
Area coordinates [L1,L2,L3] (3 * Geometrical property of the triangle (see ke_dkt) (3 * Geometrical property of the triangle (3 * Shape functions and derivatives (9 * in point L
1) 1) 1) 1)
These shape functions are used to determine the mass matrix of a triangular plate element. See also KE_DKT.
accel shell4 ACCEL_SHELL4
Compute the distributed loads for shell4 elements due to an
52
Functions — By category
acceleration. DLoads = computes In order opposite
accel_shell4(Accelxyz,Elements,Sections,Materials,Options) the distributed loads for shell4 elements due to an acceleration. to simulate gravity, accelerate the structure in the direction to gravity.
Accelxyz Elements Sections Materials Options DLoads
Acceleration Element definitions Section definitions Material definitions Element options Distributed loads
[Ax Ay Az] (1 * 3) [EltID TypID SecID MatID n1 n2 ...] [SecID SecProp1 SecProp2 ...] [MatID MatProp1 MatProp2 ... ] {Option1 Option2 ...} [EltID n1globalX n1globalY n1globalZ ...]
See also ACCEL, ACCEL_TRUSS.
loads shell4 LOADS_SHELL4
Equivalent nodal forces for a shell4 element in the GCS.
F = loads_shell4(DLoad,Node) computes the equivalent nodal forces of a distributed load (in the global coordinate system). DLoad Node F
Distributed loads in corner Nodes Node definitions Load vector (24 * 1)
[n1globalX; n1globalY; n1globalZ; ...] (12 * 1) [x y z] (4 * 3)
See also LOADSLCS_BEAM, ELEMLOADS, LOADS_TRUSS.
loadslcs shell4 LOADSLCS_SHELL4
Equivalent nodal forces for a shell4 element in the LCS.
F = loadslcs_shell4(DLoadLCS,Node) computes the equivalent nodal forces of a distributed load (in the local coordinate system). DLoadLCS Node
Distributed loads in corner Nodes Node definitions
[n1localX; n1localY; n1localZ; ...] (12 * 1) [x y z] (4 * 3)
Shell4 functions
FLCS
53
Load vector
(24 * 1)
See also LOADSLCS_BEAM, ELEMLOADS, LOADS_TRUSS.
pressure shell4 PRESSURE_SHELL4
Equivalent nodal forces for a shell4 element in the GCS due to a pressure normal to the element surface.
F = pressure_shell4(Pressure,Node) computes the equivalent nodal forces of a pressure load normal to the elements surface. Pressure
Distributed pressure in corner Nodes
Node F
Node definitions Load vector (24 * 1)
See also
[p1lobalZ; ...] (4 * 1) [x y z] (4 * 3)
ELEMPRESSURE, PRESSURE_SHELL8.
se shell4 SE_SHELL4
Compute the element stresses for a shell4 element.
[SeGCS,SeLCS,vLCS] = se_shell4(Node,Section,Material,UeGCS,Options,GCS) [SeGCS,SeLCS] = se_shell4(Node,Section,Material,UeGCS,Options,GCS) SeGCS = se_shell4(Node,Section,Material,UeGCS,Options,GCS) computes the element stresses in the global and the local coordinate system for the shell4 element. Node
Node definitions [x y z] (4 * 3) Nodes should have the following order: 4---------3 | | | | | | 1---------2
Section Material UeGCS Options
Section definition [h] Material definition [E nu rho] Displacements (24 * nTimeSteps) Element options {Option1 Option2 ...}
54
Functions — By category
GCS SeGCS
SeLCS
vLCS
Global coordinate system in which stresses are returned ’cart’|’cyl’|’sph’ Element stresses in GCS in corner nodes IJKL and at top/mid/bot of shell (72 * nTimeSteps) 72 = 6 stress comp. * 4 nodes * 3 locations [sxx syy szz sxy syz sxz] Element stresses in LCS in corner nodes IJKL and at top/mid/bot of shell (72 * nTimeSteps) [sxx syy szz sxy syz sxz] Unit vectors of LCS (1 * 9)
See also ELEMSTRESS, SE_SHELL8.
selcs shell4 SELCS_SHELL4
Compute the element stresses for a shell4 element.
[SeLCS] = selcs_shell4(Node,Section,Material,UeGCS,Options) computes the element stresses in the global and the local coordinate system for the shell4 element. Node h E nu rho UeLCS Options SeLCS
Node definitions [x y z] (4 * 3) Shell thickness Young’s modulus Poisson coefficient Mass density Displacements (24 * nTimeSteps) Element options {Option1 Option2 ...} Element stresses in LCS in corner nodes IJKL and at top/mid/bot of shell (72 * nTimeSteps) [sxx syy szz sxy syz sxz]
See also ELEMSTRESS, SE_SHELL4.
coord shell4 COORD_SHELL4
Coordinates of the shell elements for plotting.
[X,Y,Z] = coord_shell4(Nodes,NodeNum) returns the coordinates of the shell4 elements for plotting. Nodes
Node definitions
[NodID x y z] (nNodes * 4)
Shell4 functions
NodeNum X Y Z
55
Node numbers X coordinates Y coordinates Z coordinates
[NodID1 NodID2 NodID3 NodID4] (nElem * 4) (4 * nElem) (4 * nElem) (4 * nElem)
See also COORD_TRUSS, PLOTELEM.
disp shell4 DISP_SHELL4
Matrices to compute the displacements of the deformed shell4.
[Ax,Ay,Az,B] = disp_shell4(Nodes,Elements,DOF,U) returns the matrices to compute the displacements of the deformed shell4. The coordinates of the nodes of the shell4 element are computed using X=Ax*U+B(:,1); Y=Ay*U+B(:,2) and Z=Az*U+B(:,3). Nodes Elements DOF Ax Ay Az B
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Degrees of freedom (nDOF * 1) Matrix to compute the x-coordinates of the deformations Matrix to compute the y-coordinates of the deformations Matrix to compute the z-coordinates of the deformations Matrix which contains the x-, y- and z-coordinates of the undeformed structure
See also DISP_TRUSS, PLOTDISP, DISP_SHELL8.
scontour shell4 SCONTOUR_SHELL4
Matrix to plot contours in a shell4 element.
Scontour = scontour_shell4(Node,Stress,Svalues) returns a matrix used for plotting contours in a shell4 element. Node Stress Svalues Scontour
Node definitions [x y z] (4 * 3) Stress in nodes (4 * 1) Values of contours (nContour * 1) Coordinates of contours in GCS
See also PLOTSTRESSCONTOUR, PLOTSHELLFCONTOUR.
56
Functions — By category
patch shell4 PATCH_SHELL4
Patch information of the shell4 elements for plotting.
[pxyz,pind,pvalue] = patch_shell4(Nodes,NodeNum,Values) returns matrices to plot patches of shell4 elements. Nodes NodeNum Values pxyz pind pvalue
Node definitions [NodID x y z] Node numbers [NodID1 NodID2 NodID3 NodID4] (nElem * 4) Values assigned to nodes used for coloring (nElem * 4) Coordinates of Nodes (4*nElem * 3) indices of Nodes (nElem * 4) Values arranged per Node (4*nElem * 1)
See also PLOTSTRESSCONTOURF, PLOTSHELLFCONTOURF.
grid shell4 GRID_SHELL4
Grid in natural coordinates for mapped meshing.
[s,t,NodeNum,Elements] = grid_shell4(m,n,Type,Section,Material) returns matrices of a grid in the natural coordinate system, which can be used for mapped meshing. s t NodeNum Elements Type Section Material
s-coordinate of nodes (1 * nNodes) t-coordinate of nodes (1 * nNodes) Node numbers order on grid ((m+1) * (n+1)) Node numbers are saved per element here (nElem * 8) Type ID of meshed elements Section ID of meshed elements Material ID of meshed elements
See also MAKEMESH, GRID_SHELL8.
Shell8 functions
1.8
57
Shell8 functions
dof_shell8 ke_shell8 sh_qs8 b_shell8 se_shell8 accel_shell8 loads_shell8 pressure_shell8 coord_shell8 disp_shell8 scontour_shell8 patch_shell8 grid_shell8
Element degrees of freedom for a shell8 element. shell8 element stiffness and mass matrix in global coordinate system. Shape functions for a quadrilateral serendipity element with 8 nodes. B matrix for a shell8 element in global coordinate system. Compute the element stresses for a shell8 element. Compute the distributed loads for a shell due to an acceleration. Equivalent nodal forces for a shell8 element in the GCS. Equivalent nodal forces for a shell8 element in the GCS due to a pressure. Coordinates of the shell8 elements for plotting. Matrices to compute the displacements of the deformed shell. Matrix to plot contours in a shell8 element. Patch information of the shell8 elements for plotting. Grid in natural coordinates for mapped meshing.
p.58 p.58 p.59 p.59 p.60 p.61 p.61 p.61 p.62 p.62 p.63 p.63 p.63
58
Functions — By category
dof shell8 DOF_SHELL8
Element degrees of freedom for a shell8 element.
dof = dof_shell8(NodeNum) builds the vector with the labels of the degrees of freedom for which stiffness is present in the shell8 element. NodeNum Node definitions dof Degrees of freedom
[NodID1 NodID2 ... NodIDn]
(1 * 8) (48 * 1)
See also GETDOF.
ke shell8 KE_SHELL8
shell element stiffness and mass matrix in global coordinate system.
[Ke,Me] = ke_shell8(Node,Section,Material,Options) Ke = ke_shell8(Node,Section,Material,Options) returns the element stiffness and mass matrix in the global coordinate system for an eight node shell element. Node
Node definitions [x y z] (8 * 3) Nodes should have the following order: 4----7----3 | | 8 6 | | 1----5----2
Section
Section definition [h] or [h1 h2 h3 h4] (uniform thickness or defined in corner nodes(1,2,3,4)) Material definition [E nu rho] or [Exx Eyy nuxy muxy muyz muzx theta rho Element options struct. Fields: -LCSType: determine the reference local element coordinate system. Values: ’element’ (default) or ’global’ -MatType: ’isotropic’ (default) or ’orthotropic’ -Offset: nodal offset from shell midplane. Values: ’top’, ’mid’ (default), ’bot’ or numerical value -RotaryInertia: 0 (default) or 1 Element stiffness matrix (48 * 48) Element mass matrix (48 * 48)
Material Options
Ke Me
This element is based on chapter 15 of The Finite Element Method: for Solid and Structural Mechanics,
Shell8 functions
59
Zienkiewicz (2005). See also KE_BEAM, KE_BEAM, ASMKM, KE_TRUSS. KE_TRUSS.
sh qs8 SH_Q SH_QS8 S8
Shap Shape e func functi tion ons s for for an 8 node node quadr quadril ilat ater eral al ser seren endi dipi pity ty ele eleme ment nt. .
[Ni,dN_dx [Ni,dN_dxi,dN i,dN_deta _deta] ] = sh_qs8(xi sh_qs8(xi,eta ,eta) ) returns returns the shape functions functions and its deriva derivativ tives es with with respec respect t to the natura natural l coordi coordinat nates es in the point point (xi,et (xi,eta). a). xi et a Ni dN _ dx i dN _ de t a
n at u ra l co o rd i na t e n at u ra l co o rd i na t e s ha p e f u nc t io n s i n p oi n t ( x i, e ta ) de d erivative of Ni to xi derivative of Ni to eta
(s c al a r) (s c al a r) (8 * 1) (8 * 1) (8 * 1)
see also KE_SHELL8 KE_SHELL8. .
b shell8 B_SH B_SHEL ELL8 L8
b matr matrix ix for a shel shell8 l8 eleme element nt in glob global al coord coordin inat ate e system system. .
[Bg,J] = b_shell8(Ni,dN_dxi,dN_deta,z b_shell8(Ni,dN_dxi,dN_deta,zetar,Node, etar,Node,h,v1i,v2i,v h,v1i,v2i,v3i) 3i) returns the element element b matrix matrix in the global global coordina coordinate te system system and the Jacobian Jacobian of the parame parametri tric c transf transform ormati ation. on. Both Both are evalua evaluated ted in the natura natural l coordi coordinat nates es (xi,et (xi,eta a and zetar) zetar) which which were were used used to calcul calculate ate Ni,dN_dxi Ni,dN_dxi,dN_ ,dN_deta deta and zetar. zetar. No d e
Ni dN _ dx i dN _ de t a h v(1, v(1,2, 2,3) 3)i i d
Node definitions [x y z] (8 Nodes Nodes should should have have the follow following ing order: order: 4----7----3 | | 8 6 | | 1----5----2 Shape functions for quadratic serendipity element (8 first derivatives of shape functions Ni (8 first derivatives of shape functions Ni (8 scalar or vector containing thickness sc a la r o r ( 8 comp compon onen ents ts of the local local coord coordin inat ate e syste system m in node node i (8 Nodal offset from shell mid plane sc a la r ( d ef a ul t
* 3)
* * * * * =
1) 1) 1) 1) 3) 0)
60
Functions — By category
Bg J
b matrix of shell8 element Jacobian of the parametric transformation
(6 * 48) (3 * 3)
See also SE_SHELL8 SE_SHELL8, , KE_SHELL8 KE_SHELL8. .
se shell8 SE_S SE_SHE HELL LL8 8
Comp Comput ute e the the elem elemen ent t stre stress sses es for a shel shell8 l8 elemen element. t.
[SeGCS,SeLCS,vLCS] [SeGCS,SeLCS,vLCS] = se_shell8(Node,Section,Materi se_shell8(Node,Section,Material,UeGCS,O al,UeGCS,Options,gc ptions,gcs) s) [SeGCS [SeGCS,Se ,SeLCS LCS] ] = se_she se_shell8 ll8(No (Node, de,Sec Sectio tion, n,Mat Materi erial, al,UeG UeGCS, CS,Opt Option ions,g s,gcs cs) ) SeG SeGCS = se_s se_sh hell ell8(N 8(Node ode,Se ,Secti ction, on,Mat Materi erial,U l,UeGCS eGCS,O ,Opt pti ions ons,gc ,gcs) comput computes es the elemen element t stress stresses es in the global global and the local local coordi coordinat nate e system system for the shell8 shell8 elemen element. t. No d e
Node definitions [x y z] (8 * 3) Nodes Nodes should should have have the follow following ing order: order: 4----7----3 | | 8 6 | | 1----5----2
Se c ti o n
Section definition [h ] o r [ h 1 h 2 h 3 h 4 ] (uniform (uniform thickness thickness or defined defined in corner corner nodes(1,2 nodes(1,2,3,4 ,3,4)) )) Mat Materi erial def defini initio tion [E nu rho] rho] or [Ex [Exx Eyy nux nuxy muxy muxy muy muyz muzx uzx thet theta a rho rho Dis Displa placem cement ents (48 (48 * nTim nTime eSte Steps) ps) Elem Elemen ent t opti option ons s stru struct ct. . Fiel Fields ds: : -LCSType: -LCSType: determine determine the reference reference local element element coordinate coordinate system. system. Values: Values: ’element’ ’element’ (default) (default) or ’global’ ’global’ -MatType: -MatType: ’isotropic ’isotropic’ ’ (default) (default) or ’orthotrop ’orthotropic’ ic’ -Offse -Offset: t: nodal nodal offset offset from from shell shell midpla midplane. ne. Values Values: : ’top’, ’top’, ’mid’ ’mid’ (defau (default) lt), , ’bot’ ’bot’ or numeri numerical cal value value Glo Global bal coor coord dina inate syst ystem in whi which str stress esses are ret return urned ’cart’|’cyl’|’sph’ Ele Elemen ment str stre esse sses in in GC GCS in in cor corn ner nod nodes IJKL JKL an and at top/mi top/mid/b d/bot ot of shell shell (72 * nTimeS nTimeStep teps) s) 72 = 6 stress comp. * 4 nodes * 3 locations [sxx [sxx syy syy szz szz sxy sxy syz syz sxz] sxz] Ele Elemen ment str stre esse sses in in LC LCS in in cor corn ner nod nodes IJKL JKL an and at top/mi top/mid/b d/bot ot of shell shell (72 * nTimeS nTimeStep teps) s) [sxx [sxx syy syy szz szz sxy sxy syz syz sxz] sxz] Unit vectors of LCS (3 * 3)
Mate ateria rial UeGC eGCS Opti Option ons s
GCS SeGC eGCS
SeLC eLCS
vL C S
See also ELEMSTRES ELEMSTRESS, S, SE_SHELL4. SE_SHELL4.
Shell8 functions
61
accel shell8 ACCEL_ ACCEL_SHE SHELL8 LL8 DLoads = comput computes es In order order opposite opposite
Comput Compute e the distri distribut buted ed loads loads for shell8 shell8 elemen elements ts due to an accelera acceleratio tion. n.
accel_shell8(Accelxyz,Elemen accel_shell8(Accelxyz,Elements,Section ts,Sections,Materials s,Materials,Options) ,Options) the distri distribut buted ed loads loads for shell8 shell8 elemen elements ts due to an accele accelerat ration ion. . to simula simulate te gravit gravity, y, accele accelerat rate e the struct structure ure in the direct direction ion to gravity. gravity.
Accelxyz Elem Elemen ents ts Sect Sectio ions ns Mate Materi rial als s Opti Option ons s DLoa Loads
A c ce l er a ti o n [A x A y A z ] ( 1 * 3) Elem Elemen ent t defi defini niti tion ons s [Elt [EltID ID TypI TypID D SecI SecID D MatI MatID D n1 n1 n2 n2 ... ...] ] Sect Sectio ion n defi defini niti tion ons s [Sec [SecID ID SecP SecPro rop1 p1 SecP SecPro rop2 p2 ...] ...] Mate Materi rial al defi defini niti tion ons s [Mat [MatID ID MatP MatPro rop1 p1 MatP MatPro rop2 p2 ... ... ] Elem Elemen ent t opti option ons s stru struct ct. . Fiel Fields ds: : -MatType: -MatType: ’isotropic ’isotropic’ ’ (default) (default) or ’orthotrop ’orthotropic’ ic’ Dis Distri tribut buted loa loads [Elt EltID n1g n1glob lobalX alX n1g n1glob lobalY alY n1gl n1glo obal balZ ...] ...]
See also ACCEL, ACCEL, ACCEL_TRU ACCEL_TRUSS. SS.
loads shell8 LOAD LOADS_ S_SH SHEL ELL8 L8
Equi Equiva vale lent nt nodal nodal forces forces for a shel shell8 l8 eleme element nt in the the GCS. GCS.
F = loads_she loads_shell8(D ll8(DLoad Load,Node ,Node) ) comput computes es the equiva equivalen lent t nodal nodal forces forces of a distri distribut buted ed load load (in the global global coordinate coordinate system). system). DLoa Load No d e F
Dis Distri tribut buted loa loads i n c o rn e r N od e s Node definitions Load vector (4 (48 * 1)
[n1g [n1gl loba obalX; lX; n1gl n1glo obal balY; n1gl 1globa obalZ; lZ; ... ...] (1 2 * 1) [x y z] (8 * 3)
See also LOADSLCS_ LOADSLCS_BEAM BEAM, , ELEMLOADS, ELEMLOADS, LOADS_TRU LOADS_TRUSS. SS.
pressure shell8 PRES PRESSU SURE RE_S _SHE HELL LL8 8
Equi Equiva vale lent nt nodal nodal forc forces es for for a shel shell8 l8 elemen element t in the the GCS GCS due to a pressu pressure re normal normal to the element element surfac surface. e.
F = pressure_ pressure_shell shell(Pre (Pressure ssure,Nod ,Node) e) comput computes es the equiva equivalen lent t nodal nodal forces forces of a pressu pressure re load load normal normal to the elements elements surface. surface.
62
Functions — By category
Pressure
Distributed loads in corner Nodes
[p1localZ; ...] (4 * 1)
Node F See also
Node definitions Load vector (48 * 1)
[x y z] (8 * 3)
ELEMPRESSURE, PRESSURE_SHELL4.
coord shell8 COORD_SHELL8
Coordinates of the shell8 elements for plotting.
[X,Y,Z] = coord_shell8(Nodes,NodeNum) returns the coordinates of the shell8 elements for plotting. Nodes NodeNum X Y Z
Node definitions [NodID x y z] (nNodes * 4) Node numbers [NodID1 NodID2 NodID3 NodID4] (nElem * 8) X coordinates (20 * nElem) Y coordinates (20 * nElem) Z coordinates (20 * nElem)
See also COORD_TRUSS, PLOTELEM.
disp shell8 DISP_SHELL8
Matrices to compute the displacements of the deformed shell.
[Ax,Ay,Az,B] = disp_shell8(Nodes,Elements,DOF,U) returns the matrices to compute the displacements of the deformed shell. The coordinates of the nodes of the shell8 element are computed using X=Ax*U+B(:,1); Y=Ay*U+B(:,2) and Z=Az*U+B(:,3). Nodes Elements DOF Ax Ay Az B
Node definitions [NodID x y z] Element definitions [EltID TypID SecID MatID n1 n2 ...] Degrees of freedom (nDOF * 1) Matrix to compute the x-coordinates of the deformations Matrix to compute the y-coordinates of the deformations Matrix to compute the z-coordinates of the deformations Matrix which contains the x-, y- and z-coordinates of the undeformed structure
Shell8 functions
63
See also DISP_TRUSS, PLOTDISP, DISP_SHELL4.
scontour shell8 SCONTOUR_SHELL8
Matrix to plot contours in a shell8 element.
Scontour = scontour_shell8(Node,Stress,Svalues) returns a matrix used for plotting contours in a shell4 element. Node Stress Svalues Scontour
Node definitions [x y z] (4 * 3) Stress in nodes (4 * 1) or (8 * 1) Values of contours (nContour * 1) Coordinates of contours in GCS
See also PLOTSTRESSCONTOUR, PLOTSHELLFCONTOUR.
patch shell8 PATCH_SHELL8
Patch information of the shell8 elements for plotting.
[pxyz,pind,pvalue] = patch_shell8(Nodes,NodeNum,Values) returns matrices to plot patches of shell8 elements. Nodes NodeNum Values pxyz pind pvalue
Node definitions [NodID x y z] Node numbers [NodID1 NodID2 NodID3 NodID4] (nElem * 8) Values assigned to nodes used for coloring (nElem * 8) Coordinates of Nodes (8*nElem * 3) indices of Nodes (nElem * 8) Values arranged per Node (8*nElem * 1)
See also PLOTSTRESSCONTOURF, PLOTSHELLFCONTOURF.
grid shell8 GRID_SHELL8
Grid in natural coordinates for mapped meshing.
[s,t,NodeNum,Elements] = grid_shell8(m,n,Type,Section,Material) returns matrices of a grid in the natural coordinate system, which can
64
Functions — By category
be used for mapped meshing. s t NodeNum Elements Type Section Material
s-coordinate of nodes (1 * nNodes) t-coordinate of nodes (1 * nNodes) Node numbers order on grid ((m+1) * (n+1)) Node numbers are saved per element here (nElem * 12) Type ID of meshed elements Section ID of meshed elements Material ID of meshed elements
See also MAKEMESH, GRID_SHELL4.
Nonlinear analysis functions
1.9
65
Nonlinear analysis functions
nlsolver lincrpoint nlincrpoint nlincrpoint2
Solve static nonlinear problems. Interpolation-based linear estimate of critical p oints. Determine nonlinear critical points (bisection). Determine nonlinear critical points (direct method).
p.66 p.67 p.67 p.68
66
Functions — By category
nlsolver NLSOLVER
Solve static nonlinear problems
[U,K,lambda] = nlsolver(Nodes,Elements,Types,Sections,Materials,DOF,Pf); [U,K,lambda,ds,HistPar,crpoint_info] = nlsolver(Nodes,Elements,Types,... Sections,Materials,DOF,Pf,U0,Options,HistPar,Constr,P0); This function contains several algorithms for solving nonlinear problems of the form R(U)-lambda*Pf = 0 The standard algorithm uses a load control strategy with Newton-Raphson as nonlinear solver INPUT: Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } Sections Section definitions [SecID SecProp1 SecProp2 ...] Materials Material definitions [MatID MatProp1 MatProp2 ... ] DOF Degrees of freedom (nDOF * 1) Pf External load vector (nDOF * 1) Contains the final load vector (i.e. multiplier lambda varies from 0 to 1 in load control steps) U0 Initial displacement vector (nDOF * 1) Options Struct containing optional parameters. Fiels: .nloadincrem Number of load increments in load control strategy .method Solver for nonlinear system: ’NR’: Newton-Raphson (default) ’MNR’: Modified Newton-Raphson ’QN’: Quasi-Newton (BGFS) .tolerance Convergence tolerance on residual for nonlinear solvers (default: 1e-4) .arclength Use arclength control {’off’: load control (default) |’on’: arclength control} .arclength_strategy Select arclength strategy/constraint type: {’hyperplane’ (default)| ’cylinder’} .adapt_arclength Use heuristic scheme to adapt arclength during iterations {true (default)| false} .maxiter Maximum number of sub-iterations in each load or arclength substep .maxoiter Maximum number of outer iterations in arclength methods .lsmethod Add line-search method to load control strategy: {’no’ (default) | ’energy’ | ’armijo’} .initdl Initial arclength (default = 0) when set to 0 an initial arclength is estimated based on the equivalent number of load steps (nloadincrem) .ksym Assume symmetric tangent stiffness matrix {true | false (default)} .check_crpoint Check whether a critical point is passed in each step (true | false (default)) .printiter Print information during each iteration of the algorithm (true (default)| false) HistPar Element related history parameters Constr Impose linear constraint equations (see ADDCONSTR function)
Nonlinear analysis functions
67
P0 Initial load vector for solving systems: R(U)-lambda*Pf-P0 = 0 OUTPUT: U Displacement history (nDOF * nloadincrem) K Final tangent stiffness lambda Load multiplier history (1 * nloadincrem) ds Arclength history (1 * nloadincrem) HistPar Element related history parameters crpoint_info Struct containing information on critical points
lincrpoint LINCRPOINT Linear estimate of critical points based on exact geometrical stiffness in two equilibrium points. [lambda,phi] = lincrpoint(Nodes,Elements,Types,Sections,Materials,DOF,P); [lambda,phi] = lincrpoint(Nodes,Elements,Types,Sections,Materials,DOF,P,... nmode,lambda0,U0,Options); estimates the critical load factor (e.g. buckling/bifurcation load factor) by linear interpolation of the tangent stiffness matrix of two distinct equilibrium points (see p.630 Bathe and p.362 Belytschko). Input: Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } Sections Section definitions [SecID SecProp1 SecProp2 ...] Materials Material definitions [MatID MatProp1 MatProp2 ... ] DOF Degrees of freedom (nDOF * 1) P Load vector (nDOF * 1) Optional input: nmode number of modes (default: 10) lambda0 load factors of the two reference equilibrium points (default: [0 1]) U0 Initial displacements used in finding the equilibrium points (nDOF x 2) Options Options struct (passed to nlsolver function)
nlincrpoint NLINCRPOINT Determine nonlinear critical point (bisection) [lambda,phi] = nlincrpoint(Nodes,Elements,Types,Sections,Materials,DOF,P,... lambda0,U0,Options,HistPar,Constr); locate nonlinear critical point between two states using a simple
68
Functions — By category
bracketing algorithm. Input: Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } Sections Section definitions [SecID SecProp1 SecProp2 ...] Materials Material definitions [MatID MatProp1 MatProp2 ... ] DOF Degrees of freedom (nDOF * 1) P Load vector (nDOF * 1) lambda0 Load factors of the two reference equilibrium points (nDOF x 2) U0 Displacements of the reference equilibrium points (nDOF x 2) Options Options struct. Fields: .nlsolver Nonlinear solver strategy (loadcontrol (default)| arclength) .tolerance Convergence tolerance bisection algorithm (default: 1e-3) .maxiter Maximum number of iterations (default: 20) .verbose Print information during iterations (true (default)| false) HistPar Element related history parameters Constr Impose linear constraint equations (see ADDCONSTR function)
nlincrpoint2 NLINCRPOINT2 Determine nonlinear critical point (direct method) [lambda,phi] = lincrpoint2(Nodes,Elements,Types,Sections,Materials,DOF,P); [lambda,phi] = lincrpoint2(Nodes,Elements,Types,Sections,Materials,DOF,P,... nmode,lambda0,U0,phi0,Options); locates a nonlinear critical point by directly solving the extended system of equations. Not applicable to elasto-plastic problems. More information on this algorithm can be found in e.g. Wriggers (2008), Nonlinear Finite Element Methods (p.262). Input: Nodes Node definitions [NodID x y z] Elements Element definitions [EltID TypID SecID MatID n1 n2 ...] Types Element type definitions {TypID EltName Option1 ... } Sections Section definitions [SecID SecProp1 SecProp2 ...] Materials Material definitions [MatID MatProp1 MatProp2 ... ] DOF Degrees of freedom (nDOF * 1) P Load vector (nDOF * 1) lambda0 Initial guess of critical load factor (scalar) U0 Initial displacements used the equilibrium points (nDOF x 2) phi0 Initial gueess of critical mode (nDOF x 1) Options Options struct. Fields: .nphi0iters Number of iterations to estimate initial critical mode (default: 3) .tolerance Convergence tolerance on residual (default: 1e-7)
Nonlinear analysis functions
.maxiter .verbose
69
Maximum number of iterations (default: 20) Print information during iterations (true (default)| false)
Output: lambda phi
critical point load factor estimates buckling modes
(scalar) (ndof x 1)
kenl truss KENL_TRUSS
Geometric (exact) nonlinear truss element: internal force vector, element stiffness and mass matrix in global coordinate system.
[Re,Ke] = kenl_truss(Node,Section,Material,Ue,Options) returns the internal force vector, the element stiffness and mass matrix in the global coordinate system for a two node truss element (isotropic material) Node Section Material Ue Options HistPar ElemCache Re Ke Me
Node definitions [x1 y1 z1; x2 y2 z2] (2 * 3) Section definition [A] Material definition [E nu rho] Element displacement vector (6 * 1) Element options {Option1 Option2 ...} History parameters (not used) Element cache (not used) Element internal force vector (6 * 1) Element stiffness matrix (6 * 6) Element mass matrix (6 * 6)
See also KENLLCS_TRUSS, TRANS_TRUSS, NLASMKM.
kenl plane4 KENL_PLANE4
plane element internal force vector and tangent stiffness matrix in global coordinate system.
[Re,Ke] = kenl_plane4(Node,Section,Material,Ue,Options) Re = kenl_plane4(Node,Section,Material,Ue,Options) returns internal force vector and tangent stiffness matrix of the element in the global coordinate system for a four node plane element (isotropic material).
70
Functions — By category
Node
Node definitions [x y z] (4 * 3) plane4 only operates in the 2D xy-plane so z-coordinates should be equal to zero Nodes should have the following order: 4---------3 | | | | | | 1---------2
Section Material Options
Section definition [h] Material definition [E nu rho] Element options {Option1 Option2 ...} -only plane stress implemented Element internal force vector (8 * 1) Element stiffness matrix (8 * 8)
Re Ke
kenl kbeam KENL_TRUSS
Geometric (exact) nonlinear beam element: internal force vector, element stiffness and mass matrix in global coordinate system.
[Re,Ke] = kenl_beam(Node,Section,Material,Ue,Options) returns the internal force vector, the element stiffness and mass matrix in the global coordinate system for a two node beam element (isotropic material) Node Section Material Ue Options HistPar ElemCache Re Ke Me
Node definitions [x1 y1 z1; x2 y2 z2] (2 * 3) Section definition [A] Material definition [E nu rho] Element displacement vector (24 * 1) Element options {Option1 Option2 ...} History parameters (not used) Element cache (not used) Element internal force vector (36 * 1) Element stiffness matrix (36 * 36) Element mass matrix (36 * 36)
Implementation of beam element described in: Nielsen, M. & Krenk, S. Explicit free-floating beam element International Journal for Numerical Methods in Engineering, 2014, 98, 59-78
Chapter 2
Tutorial
72
2.1
Tutorial
Tutorial: static analysis
A frame is clamped at node 1 and has a hinge at node 5. An internal hinge is present at nodes 2 and 3. At node 4 a point load F = 5 kN is applied and the vertical beam on the left is loaded by a distributed load p = 2 kN/m. The frame consists of rectangular concrete beams (Young’s modulus E = 30 × 106 kN/m2 and Poisson coefficient ν = 0.2) with a height of 0.25m and a width of 0.1 m. A diagonal steel cable (Young’s modulus E = 210 × 106 kN/m2 and Poisson coefficient ν = 0.3) with a radius of 0.004 m, supports the frame.
Figure 2.1: Example 1: nodes, loads and boundary conditions. % StaBIL manual % Tutorial: static analysis % Units: m, kN % Nodes=[NodID Nodes= [1 2 3 4 5 6
X 0 0 0 4 4 1
Y 0 4 4 4 0 5
Z] 0; 0; 0; 0; 0; 0];
% reference node
% Check the node coordinates as follows: figure plotnodes(Nodes); % Element types -> {EltTypID EltName} Types= {1 beam ; 2 truss }; b=0.10;
Tutorial: static analysis
73
h=0.25; r=0.004; % Sections=[SecID A ky Sections= [1 b*h Inf 2 pi*r^2 NaN
kz Inf NaN
Ixx Iyy Izz 0 0 b*h^3/12 NaN NaN NaN
% Materials=[MatID E nu]; Materials= [ 1 30e6 0.2; 2 210e6 0.3]; % Elements=[EltID Elements= [1 2 3 4
TypID 1 1 1 2
SecID 1 1 1 2
MatID 1 1 1 2
yt h/2 NaN
yb h/2 NaN
zt b/2 NaN
zb] b/2; NaN];
% concrete % steel n1 1 3 5 1
n2 2 4 4 4
n3] 6; 6; 6; NaN];
% Check node and element definitions as follows: hold( on ); plotelem(Nodes,Elements,Types); title( Nodes and elements ); % Degrees of freedom % Assemble a column matrix containing all DOFs at which stiffness is % present in the model: DOF=getdof(Elements,Types); % Remove all DOFs equal to zero from the vector: % - 2D analysis: select only UX,UY,ROTZ % - clamp node 1 % - hinge at node 5 seldof=[0.03; 0.04; 0.05; 1.00; 5.01; 5.02]; DOF=removedof(DOF,seldof); % Assembly of stiffness matrix K K=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % Nodal loads: 5 kN horizontally on node 4. seldof=[4.01]; PLoad= [5]; % Assembly of the load vectors: P=nodalvalues(DOF,seldof,PLoad); % Distributed loads are specified in the global coordinate system % DLoads=[EltID n1globalX n1globalY n1globalZ ...] DLoads= [1 2 0 0 2 0 0]; P=P+elemloads(DLoads,Nodes,Elements,Types,DOF); % Constraint equations: Constant=Coef1*DOF1+Coef2*DOF2+ ... % Constraints=[Constant Coef1 DOF1 Coef2 DOF2 ...] Constr= [0 1 2.01 -1 3.01; 0 1 2.02 -1 3.02]; % Add constraint equations
74
[K,P]=addconstr(Constr,DOF,K,P); % Solve K * U = P U=K\P; % Plot displacements figure plotdisp(Nodes,Elements,Types,DOF,U,DLoads,Sections,Materials) % The displacements can be displayed as follows: printdisp(Nodes,DOF,U); % Compute element forces Forces=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads); % The element forces can be displayed in a orderly table: printforc(Elements,Forces); % Plot element forces figure plotforc( norm ,Nodes,Elements,Types,Forces,DLoads) title( Normal forces ) figure plotforc( sheary ,Nodes,Elements,Types,Forces,DLoads) title( Shear forces ) figure plotforc( momz ,Nodes,Elements,Types,Forces,DLoads) title( Bending moments ) % Plot stresses figure plotstress( snorm ,Nodes,Elements,Types,Sections,Forces,DLoads) title( Normal stresses due to normal forces ) figure plotstress( smomzt ,Nodes,Elements,Types,Sections,Forces,DLoads) title( Normal stresses due to bending moments around z: top ) figure plotstress( smomzb ,Nodes,Elements,Types,Sections,Forces,DLoads) title( Normal stresses due to bending moments around z: bottom ) figure plotstress( smax ,Nodes,Elements,Types,Sections,Forces,DLoads) title( Maximal normal stresses ) figure plotstress( smin ,Nodes,Elements,Types,Sections,Forces,DLoads) title( Minimal normal stresses )
Tutorial
76
Tutorial
173714
−71 24
−271 −71
24
−271
173714
(a) normal forces 2325 −2325
0
0
760
−760
−8525
0
(b) b ending moments around z: top
8525
0
(c) b ending moments around z: b ottom 173714
−71 24
2325 −2325
173714
2254 2054
−71 24
785
−2597 −2397
−736
173714
173714
8549
(c) maximal
−271
−8501
(d) minimal
Figure 2.3: Normal stresses.
−271
Tutorial: dynamic analysis
2.2
77
Tutorial: dynamic analysis
A dynamic analysis is made of the frame that was treated in section 2.1. % StaBIL manual % Tutorial: dynamic analysis: model % Units: m, N L=4; H=4; nElemCable=8; % Nodes=[NodID X Y Z] Nodes= [1 0 0 0; 2 L 0 0]; Nodes=reprow(Nodes,1:2,4,[2 0 H/4 0]); Nodes= [Nodes; 11 0 H 0]; Nodes=reprow(Nodes,11,3,[1 L/4 0 0]); Nodes= [Nodes; 15 1 5 0]; % reference node Nodes= [Nodes; 16 0 0 0]; Nodes=reprow(Nodes,16,nElemCable,[1 L/nElemCable H/nElemCable 0]); % Check the node coordinates as follows: figure plotnodes(Nodes); % Element types -> {EltTypID EltName} Types= {1 beam }; b=0.10; h=0.25; r=0.004; % Sections=[SecID A ky Sections= [1 b*h Inf 2 pi*r^2 Inf
kz Inf Inf
Ixx Iyy Izz] 0 0 b*h^3/12; 0 0 pi*r^4/4];
% Materials=[MatID E nu]; Materials= [1 30e9 0.2 2500; 2 210e9 0.3 7850];
% concrete % steel
% Elements=[EltID TypID SecID MatID n1 n2 n3] Elements= [1 1 1 1 1 3 15; 2 1 1 1 2 4 15]; Elements=reprow(Elements,1:2,3,[2 0 0 0 2 2 0]); Elements=[ Elements; 9 1 1 1 11 12 15; 10 1 1 1 12 13 15; 11 1 1 1 13 14 15; 12 1 1 1 14 10 15; 13 1 2 2 16 17 15]; Elements=reprow(Elements,13,(nElemCable-1),[1 0 0 0 1 1 0]);
78
Tutorial
% Check node and element definitions as follows: hold( on ); plotelem(Nodes,Elements,Types); title( Nodes and elements ); % Degrees of freedom DOF=getdof(Elements,Types); % Boundary conditions seldof=[0.03; 0.04; 0.05; 1.01; 1.02; 1.06; 2.01; 2.02; 16.01; 16.02]; DOF=removedof(DOF,seldof); % Assembly of stiffness matrix K [K,M]=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % DLoads=[EltID n1globalX n1globalY n1globalZ ...] DLoads=[1 2000 0 0 2000 0 0; 3 2000 0 0 2000 0 0; 5 2000 0 0 2000 0 0; 7 2000 0 0 2000 0 0]; b=elemloads(DLoads,Nodes,Elements,Types,DOF); % Spatial distribution, nodal (nDOF * 1) % Constraint equations: Constant=Coef1*DOF1+Coef2*DOF2+ ... % Constraints=[Constant Coef1 DOF1 Coef2 DOF2 ...] Constr= [0 1 9.01 -1 11.01; 0 1 9.02 -1 11.02; 0 1 10.01 -1 (16.01+nElemCable); 0 1 10.02 -1 (16.02+nElemCable)]; [K,b,M]=addconstr(Constr,DOF,K,b,M);
2.2.1
Eigenmodes
% StaBIL manual % Tutorial: dynamic analysis: eigenvalue problem % Units: m, N % Assembly of M and K tutorialdyna; % Eigenvalue problem nMode=12; [phi,omega]=eigfem(K,M,nMode); % Display eigenfrequenties disp( Lowest eigenfrequencies [Hz] ); disp(omega/2/pi); % Plot eigenmodes figure; plotdisp(Nodes,Elements,Types,DOF,phi(:,1), DispMax , off ) figure; plotdisp(Nodes,Elements,Types,DOF,phi(:,2), DispMax , off ) figure;
Tutorial: dynamic analysis
79
plotdisp(Nodes,Elements,Types,DOF,phi(:,5), DispMax , off ) figure; plotdisp(Nodes,Elements,Types,DOF,phi(:,8), DispMax , off ) figure; plotdisp(Nodes,Elements,Types,DOF,phi(:,11), DispMax , off ) figure; plotdisp(Nodes,Elements,Types,DOF,phi(:,12), DispMax , off ) % Animate eigenmodes figure; animdisp(Nodes,Elements,Types,DOF,phi(:,1)) title( Eigenmode 1 ) figure; animdisp(Nodes,Elements,Types,DOF,phi(:,2)) title( Eigenmode 2 ) figure; animdisp(Nodes,Elements,Types,DOF,phi(:,5)) title( Eigenmode 5 ) figure; animdisp(Nodes,Elements,Types,DOF,phi(:,8)) title( Eigenmode 8 ) figure; animdisp(Nodes,Elements,Types,DOF,phi(:,11)) title( Eigenmode 11 ) figure; animdisp(Nodes,Elements,Types,DOF,phi(:,12)) title( Eigenmode 12 )
2.2.2
Modal superposition: time domain: piecewise exact integration
% StaBIL manual % Tutorial: dynamic analysis: modal superposition: piecewise exact integration % Units: m, N % Assembly of M and K tutorialdyna; % Sampling parameters T=2.5; dt=0.002; N=T/dt; t=[0:N-1]*dt;
% % % %
% Eigenvalue analysis nMode=12; [phi,omega]=eigfem(K,M,nMode); xi=0.07;
% Number of modes to take into account % Calculate eigenmodes and eigenfrequencies % Constant modal damping ratio
% Excitation bm=phi. *b;
% Spatial distribution, modal (nMode * 1)
Time window Time step Number of samples Time axis
80
Tutorial
(a) Eigenmode 1
(b) Eigenmode 2
(b) Eigenmode 5
(c) Eigenmode 8
(d) Eigenmode 11
(e) Eigenmode 12
Figure 2.4: Results of example 1.
Tutorial: dynamic analysis
q=zeros(1,N); q((t>=0.50) & (t<0.60))=1; pm=bm*q;
81
% Time history (1 * N) % Time history (1 * N) % Modal excitation (nMode * N)
% Modal analysis x= msupt(omega,xi,t,pm, zoh ); % Modal displacements -> nodal displacements u=phi*x; % Nodal response (nDOF * N) % Figures figure; plot(t,x); title( Modal response (piecewise linear exact integration) ); xlabel( Time [s] ); xlim([0 4.1]) ylabel( Displacement [m kg^{0.5}] ); legend([repmat( Mode ,nMode,1) num2str([1:nMode]. )]); figure; c=selectdof(DOF,[9.01; 13.02; 17.02]); plot(t,c*u); title( Nodal response (piecewise linear exact integration) ); xlabel( Time [s] ); xlim([0 4.1]) ylabel( Displacement [m] ); legend( 9.01 , 13.02 , 17.02 ); % Movie figure; animdisp(Nodes,Elements,Types,DOF,u); % Display disp( Maximum modal response ); disp( max(abs(x),[],2)); disp( Maximum nodal response 9.01 13.02 17.02 ); disp( max(abs(c*u),[],2));
2.2.3
Modal superposition: transform to frequency domain
% StaBIL manual % Tutorial: dynamic analysis: direct method: frequency domain % Units: m, N % Assembly of M and K tutorialdyna; % Sampling parameters N=2048; dt=0.002; T=N*dt; F=N/T; df=1/T; t=[0:N-1]*dt;
% % % % % %
Number of samples Time step Period Sampling frequency Frequency resolution Time axis
82
Tutorial
f=[0:N/2-1]*df; Omega=2*pi Omega=2*pi*f; *f;
% Pos Positiv itive e fre freque quenci ncies es cor corres respon pondin ding g to FFT [Hz [Hz] ] % Ide Idem m [ra [rad/ d/s] s]
% Eigenv Eigenvalue alue analysis nMode=12; [phi,omega]=eigfem [phi,omega]=eigfem(K,M,nMode); (K,M,nMode); xi=0.07;
% Numb Number er of mo mode des s to ta take ke in into to ac acco coun unt t % Calcula Calculate te eigenmodes eigenmodes and eigenf eigenfrequen requencies cies % Con Constan stant t mod modal al dam dampin ping g rat ratio io
% Excita Excitation tion bm=phi. *b; q=zeros q=zeros(1,N); (1,N); q((t>=0.50) q((t>=0.50) & (t<0.60))=1; (t<0.60))=1; Q=fft Q=fft(q); (q); Q=Q(1:N/2); Pm=bm*Q;
% % % % % %
Spatial dis Spatial distri tribut bution ion, , mod modal al (nM (nMode ode * 1) Time Ti me hi hist stor ory y (1 * N) Time Ti me hi hist stor ory y (1 * N) Freq Fr equen uency cy co cont nten ent t (1 * N) Frequen Fre quency cy con conten tent, t, pos positi itive ve fre freq q (1 * N/2 N/2) ) Modal Mod al exc excita itatio tion n, pos positiv itive e fre freq q (nM (nMode ode * N/2 N/2) )
% Mod Modal al ana analys lysis is [X,H]= msupf msupf(omega,xi,Omega,Pm); (omega,xi,Omega,Pm); % Mod Modal al res respon ponse se, , pos positiv itive e fre freq q (nM (nMode ode * N/2 N/2) ) % FF-do dom m -> tt-do dom m X=[X, zeros(nMode,1), zeros(nMode,1), conj(X(:, conj(X(:,end end:-1:2))]; :-1:2))]; x=ifft x=ifft(X,[],2); (X,[],2); % Mo Moda dal l re resp spon onse se (n (nMo Mode de * N) % Mod Modal al dis displa placem cement ents s -> nod nodal al dis displa placem cements ents u=phi*x; % No Noda dal l re resp spon onse se (n (nDO DOF F * N) % Fig Figure ures s figure; figure ; subplot(3,2,1); subplot (3,2,1); plot(t,q, plot (t,q, .- ); xlim([0 xlim([0 4.1]) ylim([0 ylim([0 1.2]); title( title ( Excitation Excitation time histor history y ); xlabel( xlabel ( Time [s] ); ylabel( ylabel ( Force [N/m] ); subplot(3,2,2); subplot(3,2,2); plot(f, plot (f,abs abs(Q)/F, (Q)/F, .- ); title( title ( Excitation Excitation frequen frequency cy conten content t ); xlabel( xlabel ( Frequency Frequency [Hz] ); ylabel( ylabel ( Force [N/m/ [N/m/Hz] Hz] ); subplot(3,2,4); subplot(3,2,4); plot(f, plot (f,abs abs(H), (H), .- ); title( title ( Modal transf transfer er functi function on ); xlabel( xlabel ( Frequency Frequency [Hz] ); ylabel( ylabel ( Displacement Displacement [m/N] ); legend([repmat( legend ([repmat( Mode ,nMode,1) num2str([1:nMode]. num2str([1:nMode]. )]); subplot(3,2,6); subplot(3,2,6); plot(f, plot (f,abs abs(X(:,1:N/2))/F, (X(:,1:N/2))/F, .- ); title( title ( Modal respon response se ); xlabel( xlabel ( Frequency Frequency [Hz] ); ylabel( ylabel ( Displacement Displacement [m kg^{0 kg^{0.5}/ .5}/Hz] Hz] ); subplot(3,2,5); subplot (3,2,5);
Tutorial: dynamic analysis
83
plot(t,x); plot(t,x); title( title ( Modal response (calculation (calculation in f-dom f-dom) ) ); xlabel( xlabel ( Time [s] ); xlim([0 xlim([0 4.1]) ylabel( ylabel ( Displacement Displacement [m kg^{0 kg^{0.5}] .5}] ); figure; figure; plot(t,x); plot (t,x); title( title ( Modal response (calculation (calculation in f-dom f-dom) ) ); xlabel( xlabel ( Time [s] ); xlim([0 xlim([0 4.1]) ylabel( ylabel ( Displacement Displacement [m kg^{0 kg^{0.5}] .5}] ); legend([repmat( legend ([repmat( Mode ,nMode,1) num2str([1:nMode]. num2str([1:nMode]. )]); figure; figure; c=selectdof c=selectdof(DOF (DOF,[9.01 ,[9.01; ; 13.02; 17.02]); 17.02]); plot(t,c*u); plot (t,c*u); title( title ( Nodal response (calculation (calculation in f-dom f-dom) ) ); xlabel( xlabel ( Time [s] ); xlim([0 xlim([0 4.1]) ylabel( ylabel ( Displacement Displacement [m] ); legend( legend ( 9.01 , 13.02 , 17.02 ); % Mo Movi vie e figure; figure ; animdisp(Nodes,Elements,Types,DOF,u); animdisp (Nodes,Elements,Types,DOF,u); % Dis Displa play y disp( disp ( Maximum Maximum modal respon response se ); disp( disp ( max max( (abs abs(x),[],2)); (x),[],2)); disp( Maximum disp( Maximum nod nodal al res respon ponse se 9.0 9.01 1 13. 13.02 02 17. 17.02 02 ); disp( disp ( max max( (abs abs(c*u),[],2)); (c*u),[],2));
2.2.4 2.2.4
Direc Directt tim time e int integr egrati ation on
% StaBIL StaBIL manual % Tutori Tutorial: al: dynamic analys analysis: is: direct time integration: integration: trapez trapezium ium rule % Un Unit its s: m, N % Assembly of M, K and C tutorialdyna; [phi,omega]=eigfem [phi,omega]=eigfem(K,M); (K,M); % Calcul Calculate ate eigenm eigenmodes odes and eigenf eigenfrequen requencies cies xi=0.07; % Dam Dampin ping g rat ratio io nModes=length nModes=length(K)(K)-size size(Constr,1); (Constr,1); C=M. *phi(:,1:nModes)*diag *phi(:,1:nModes)*diag(2*xi*omega(1:nModes))*phi(:,1:nModes). (2*xi*omega(1:nModes))*phi(:,1:nModes). *M; % Mo Moda dal l -> fu full ll da damp mpin ing g ma matr trix ix C % Sampli Sampling ng parame parameters ters T=2.5; dt=0.002; N=T/dt; t=[0:N-1]*dt; % Excita Excitation tion
% % % %
Time win Time window dow Time Ti me st step ep Number Num ber of sam sample ples s Time Ti me ax axis is
84
Tutorial
Excitation time history
Excitation frequency content 0.1
1
0.08
] z H / m0.06 / N [ e 0.04 c r o F
] 0.8 m / N [ e 0.6 c r o F 0.4
0.02
0.2 0 0
1
2 Time [s]
3
0
4
40
50
] N0.6 / m [ t n e m0.4 e c a l p s 0.2 i D
0
0.15
0
10 −3
8
20 30 Frequency [Hz]
40
50
40
50
Modal response
x 10
] z H /
]
5 . 0
0.1
6
g k m [ t 4 n e m e c a 2 l p s i D
0.05
0
−0.05 0
20 30 Frequency [Hz]
0.8
Modal response (calculation in f−dom)
g k m [ t n e m e c a l p s i D
10
Modal transfer function
Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 6 Mode 7 Mode 8 Mode 9 Mode 10 Mode 11 Mode 12
5 . 0
0
1
2 Time [s]
3
4
0
0
10
20 30 Frequency [Hz]
Figure 2.5: Modal superposition in frequency domain. q=zeros q=zeros(1,N); (1,N); q((t>=0.50) q((t>=0.50) & (t<0.60))=1; (t<0.60))=1; p=b*q;
% Ti Time me hi hist stor ory y (1 * N) % Ti Time me hi hist stor ory y (1 * N) % No Noda dal l ex exci cita tati tion on (n (nDO DOF F * N)
% Dir Direct ect time int integr egrati ation on - tra trapez pezium ium rule alpha=1/4; delta=1/2; theta=1; u=wilson u=wilson(M,C (M,C,K,dt ,K,dt,p,[a ,p,[alpha lpha delta theta]); theta]);
Tutorial: dynamic analysis
85
0.1
0.1 Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 6 Mode 7 Mode 8 Mode 9 Mode 10 Mode 11 Mode 12
0.08 ]
5 . 0
0.06
g k m [ t n e m e c a l p s i D
0.04 0.02 0
−0.02 −0.04
Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 6 Mode 7 Mode 8 Mode 9 Mode 10 Mode 11 Mode 12
0.08 ]
5 . 0
0.06
g k m [ t n e m e c a l p s i D
0.04 0.02 0
−0.02
0
0.5
1
1.5
2 2.5 Time [s]
3
3.5
(a) Modal superposition in time domain
4
−0.04 0
0.5
1
1.5
2 2.5 Time [s]
(b) Modal superposition in frequency domain
Figure 2.6: Modal response: comparison.
% Figures figure; c=selectdof(DOF,[9.01; 13.02; 17.02]); plot(t,c*u); title([ Nodal response (direct time integration) ]); xlabel( Time [s] ); xlim([0 4.1]) ylabel( Nodal displacements [m] ); legend( 9.01 , 13.02 , 17.02 ); % Movie figure; animdisp(Nodes,Elements,Types,DOF,u); % Display disp( Maximum nodal response 9.01 13.02 17.02 ); disp( max(abs(c*u),[],2));
2.2.5
Direct solution in the frequency domain
% StaBIL manual % Tutorial: dynamic analysis: modal superposition: transform to f-dom % Units: m, N % Assembly of M, K and C tutorialdyna; [phi,omega]=eigfem(K,M); % Calculate eigenmodes and eigenfrequencies xi=0.07; % Damping ratio nModes=length(K)-size(Constr,1); C=M. *phi(:,1:nModes)*diag(2*xi*omega(1:nModes))*phi(:,1:nModes). *M; % Modal -> full damping matrix C % Sampling parameters
3
3.5
4
86
Tutorial
N=2048; dt=0.002; T=N*dt; F=N/T; df=1/T; t=[0:N-1]*dt; f=[0:N/2-1]*df; Omega=2*pi*f; % Excitation q=zeros(1,N); q((t>=0.50) & (t<0.60))=1; Q=fft(q); Q=Q(1:N/2); Pd=b*Q;
% % % % % % % %
Number of samples Time step Period Sampling frequency Frequency resolution Time axis Positive frequencies corresponding to FFT [Hz] Idem [rad/s]
% Time history (1 * N) % Time history (1 * N) % Frequency content (1 * N) % Frequency content, positive freq (1 * N/2) % Nodal excitation, positive freq (nDOF * N/2)
% Solution for each frequency Ud=zeros(size(Pd)); for k=1:N/2 Kd=-Omega(k)^2*M+Omega(k)*i*C+K; Ud(:,k)=Kd\Pd(:,k); end % F-dom -> t-dom Ud=[Ud, zeros(length(K),1), conj(Ud(:,end:-1:2))]; u=ifft(Ud,[],2); % Nodal response (nDOF * N) % Figures figure; c=selectdof(DOF,[9.01; 13.02; 17.02]); plot(t,c*u); title( Nodal response (direct method in f-dom) ); xlabel( Time [s] ); xlim([0 4.1]) ylabel( Displacement [m] ); legend( 9.01 , 13.02 , 17.02 ); % Movie figure; animdisp(Nodes,Elements,Types,DOF,u); % Display disp( Maximum nodal response 9.01 13.02 17.02 ); disp( max(abs(c*u),[],2));
2.2.6
Comparison
Computational cost in seconds Modal superposition Direct method
Time domain 0.167 0.213
Freq domain 0.147 0.789
Tutorial: dynamic analysis
87
−3
5
−3
x 10
5 9.01 13.02 17.02
4
] m [ t n e m e c a l p s i D
3
1
3 2 1
0
0
−1
−1
−2
0
0.5
1
1.5
2 2.5 Time [s]
3
3.5
4
(a) Modal superposition in time domain
−2 0
5 9.01 13.02 17.02
3
1
2 2.5 Time [s]
3
3.5
4
1
−1
−1
1
1.5
(b) Direct time integration
2 2.5 Time [s]
3
3.5
4
9.01 13.02 17.02
2
0
0.5
x 10
3
0
0
1.5
4
] m [ t n e m e c a l p s i D
2
−2
1
−3
x 10
4 ] m [ s t n e m e c a l p s i d l a d o N
0.5
(b) Modal superposition in frequency domain
−3
5
9.01 13.02 17.02
4
] m [ t n e m e c a l p s i D
2
x 10
−2 0
0.5
1
1.5
2 2.5 Time [s]
(c) Direct method in frequency domain
Figure 2.7: Nodal response: comparison.
3
3.5
4
88
2.3
Tutorial
Tutorial: shell analysis
A barrel vault roof subjected to its self weight is analysed. The curved edges are simply supported and the straight edges are free. Due to symmetry only a quarter of the roof is modelled and symmetry boundary conditions are applied.
Edge 2: simply supported
Edge 1: symmetry conditions
Edge 3: free boundary
zy x
Edge 4: symmetry conditions
Figure 2.8: 8 × 8 mesh of a quarter of a cylindrical roof. % A BARREL VAULT ROOF SUBJECTED TO ITS SELF WEIGHT % Reference: COOK, CONCEPTS AND APPL OF F.E.A., 2ND ED., 1981, PP. 284-287. %% Parameters R=25; L=50; t=0.25; theta = 40*pi/180; E = 4.32*10^8; nu = 0; rho = 36.7347; N = 8;
% % % % % % % %
Radius of cylindrical roof Length of cylindrical roof Thickness of roof Angle of cylindrical roof Youngs modulus Poisson coefficient Density Number of elements
%% Mesh % Lines Line1 = Line2 = Line3 =
= [x1 y1 z1;x2 y2 z2;...] [0 0 0;0 0 L/2]; [R*sin(theta*(0:0.1:1). ) R*cos(theta*(0:0.1:1). )-R L*ones(11,1)/2]; [R*sin(theta) R*cos(theta)-R L/2; R*sin(theta) R*cos(theta)-R 0];
Tutorial: shell analysis
Line4 = [R*sin(theta*(1:-0.1:0). ) R*cos(theta*(1:-0.1:0). )-R zeros(11,1)]; % Specify element type for mesh Materials = [1 E nu rho]; Sections = [1 t]; Types = {1 shell8 }; % Mesh the area between lines 1,2,3,4 with N * N elements of type shell8, % section number 1 and material number 1. [Nodes,Elements,Edge1,Edge2,Edge3,Edge4] = ... makemesh(Line1,Line2,Line3,Line4,N,N,Types(1,:),1,1); % Check mesh: figure; plotnodes(Nodes, numbering , off ); hold( on ) plotelem(Nodes,Elements,Types, numbering , off ); title( Nodes and elements ); %% Assemble stiffness matrix % Select all dof: DOF = getdof(Elements,Types); % Apply boundary conditions: % - Line1 and Line4: symmetry condition % - Line2: simply supported sdof = [Edge1+0.01;Edge1+0.06;Edge1+0.05;Edge2+0.02;Edge2+0.01; Edge2+0.06;Edge4+0.03;Edge4+0.04;Edge4+0.05]; DOF = removedof(DOF,sdof); % Assemble K: K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF); %% Solution % Apply gravitational acceleration and determine equivalent nodal forces: DLoads=accel([0 9.8 0],Elements,Types,Sections,Materials); P=elemloads(DLoads,Nodes,Elements,Types,DOF);
% Solve K * U = P: U = K\P; % Plot displacements: figure; plotdisp(Nodes,Elements,Types,DOF,U) title( Displacements ) % Check target displacement: TP1 = selectdof(DOF,intersect(Edge3,Edge4)+0.02); Utp1 = TP1*U ratio_u = -Utp1/0.3016 %% Stress
89
90
% Determine element stress in global and local(element) coordinate system: [SeGCS,SeLCS,vLCS]=elemstress(Nodes,Elements,Types,Sections,Materials,DOF,U); % print stress: printstress(Elements,SeGCS) % plot stress contour: figure; plotstresscontour( sx ,Nodes,Elements,Types,SeGCS, location , bot ) title( sx in gcs (element solution) ) % plot filled contours: figure; plotstresscontourf( sx ,Nodes,Elements,Types,SeGCS, location , bot ) title( sx in gcs (element solution) ) figure; plotstresscontour( sx ,Nodes,Elements,Types,SeLCS, location , bot ) title( sx in lcs ) % plot lcs for shell elements figure; plotlcs(Nodes,Elements,Types,vLCS) title( local coordinate system ) % Calculate nodal solution % SnGCS: stress arranged per element % SnGCS2: stress arranged per node [SnGCS,SnGCS2]=nodalstress(Nodes,Elements,Types,SeGCS); figure; plotstresscontour( sx ,Nodes,Elements,Types,SnGCS, location , bot ); title( sx in gcs (nodal solution) ) % Stress ratios: ratio_sz = SnGCS2(intersect(Edge3,Edge4),16)/358420 ratio_st = SnGCS2(intersect(Edge4,Edge1),14)/(-213400) %% Shell forces % Shell forces (element solution): [FeLCS]=elemshellf(Elements,Sections,SeLCS); figure; plotshellfcontour( my ,Nodes,Elements,Types,FeLCS) title( my (element solution) ) % Nodal solution: [FnLCS,FnLCS2]=nodalshellf(Nodes,Elements,Types,FeLCS); figure; plotshellfcontour( my ,Nodes,Elements,Types,FnLCS) title( my (nodal solution) ) %% Principal stress % Principal stresses: [Spr,Vpr]=principalstress(Elements,SnGCS);
Tutorial
Tutorial: shell analysis
figure; plotstresscontour( s1 ,Nodes,Elements,Types,Spr, location , bot ); title( s1 ) % plot principal stresses: figure; plotprincstress(Nodes,Elements,Types,Spr,Vpr) title( principal stress directions )
91
92
Tutorial
zy x
(a) Local coordinate systems
(b) Displacements
5
0.0599
5
x 10
3.5841
−0.1401
3.2082
−0.3401
2.8323
−0.5401
2.4564
−0.7401
2.0804
−0.9401
1.7045
−1.1401
1.3286
−1.3401
0.9527 zy
−1.5401 x −1.7401 −1.9401
σxx
zy
0.5768 x
MIN
0.2008
MIN
−0.1751
MAX
−2.1401
(c)
x 10
−0.551
at bottom of shell in GCS (element solution)
(d)
σxx
MAX
at bottom of shell in LCS (element solution)
2079.965 1881.9048 1683.8446 1485.7844 1287.7242 1089.6639 891.6037 693.5435 zy
495.4833 x 297.4231 99.3629
zy
MAX
x
MIN
−98.6973
(e)
my
in LCS (nodal solution)
(f ) Principal stresses at top of shell
Figure 2.9: Cylindrical roof: results.
Tutorial: nonlinear analysis
2.4
93
Tutorial: nonlinear analysis
This tutorial performs a geometric nonlinear analysis of an arc loaded by a point force. Several methods are used to compute the critical buckling load. clear %% GEOMETRIC NONLINEAR ARC % (Example based on p268 of Rin = 100; Rout = 103; alpha = 60*pi/180; E = 40000; nu = 0.2; m = 6; n = 150; Types = {1 plane4 }; Sections = [1 1]; Materials = [1 E nu];
Nonlinear Finite Element Methods, Wriggers (2008)) % Inner radius arc % Outer radius arc % Arc angle % Young s modulus % Poisson s coefficient % number of elements in radial direction % number of elements in tangential direction % {EltTypID EltName} % [SecID t] % [MatID E nu]
%% MESH npoints = 100; theta = linspace(-alpha/2+pi/2,alpha/2+pi/2,npoints). ; [Xin,Yin] = pol2cart(flipud(theta),Rin*ones(npoints,1)); [Xout,Yout] = pol2cart(theta,Rout*ones(npoints,1)); Yin = Yin-Rin*cos(alpha/2); Yout = Yout-Rin*cos(alpha/2); Z = zeros(npoints,1); Line1 = [Xin(end) Yin(end) 0;Xout(1) Yout(1) 0]; Line2 = [Xout Yout Z]; Line3 = [Xout(end) Yout(end) 0;Xin(1) Yin(1) 0]; Line4 = [Xin Yin Z]; [Nodes,Elements,Edge1,Edge2,Edge3,Edge4] = makemesh(Line1,Line2,Line3,Line4,n,m,Types(1,:),1,1); figure; plotelem(Nodes,Elements,Types, numbering , off ); title( Elements ); %% DOF DOF = getdof(Elements,Types); sdof = [Edge1((end-1)/2+1);Edge3((end-1)/2+1)]; sdof = [sdof+0.01;sdof+0.02]; DOF = removedof(DOF,sdof); %% LOAD pdof = Edge2((end-1)/2+1)+0.02; Lp = selectdof(DOF,pdof); pmax = 1000; P = -pmax*Lp. ; %% LINEAR ANALYSIS K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF); U = K\P; %% Newton Raphson Options.method = nr ; Options.nloadincrem = 50;
% Nonlinear solver: Newton-Raphson % Number of load increments
Tutorial: nonlinear analysis
95
hold on; plot(-Lp*Ual,lambda_al*pmax); plot(-Lp*Ual(:,alstability.nonpositivepivots>0),lambda_al(alstability.nonpositivepivots>0)*pmax, x- ) xlabel( U ); ylabel( P );axis tight; ax = axis; plot([ax(1),ax(2)],[lambda(1) lambda(1)], k-- ); plot(-Lp*Unlcr,lambdanlcr*pmax, rx ); plot(-Lp*Ual2,lambda_al2*pmax, r-- ) legend( Stable path , Unstable path , Linear critical load ,... Nonlinear critical load , Bifurcated path ); title( Load-displacement curve ); hold off; %% ANIMATE LOAD DISPLACEMENT CURVE nstep = size(Ual,2); figure for istep = [(1:10:nstep),nstep] clf subplot(2,1,1) plotdisp(Nodes,Elements,Types,DOF,Ual(:,istep), dispscal ,1, dispmax , off ); subplot(2,1,2) plot(-Lp*Ual(:,1:istep),pmax*lambda_al(:,1:istep), .- ); xlabel( U ); ylabel( P ); axis(ax); pause(0.1) end
96
Tutorial
1000
Stable path Unstable path Linear critical load Nonlinear critical load Bifurcated path
800 600 400 200 P
y z x
0 −200 −400 −600 −800 0
(a) Mesh
Maximal displacement: 8.76453
(c) Linear small displacements
5
10
15 U
(b) Nonlinear load-displacement curve
Maximal displacement: 30.7559
(d) Nonlinear displacements
Figure 2.10: Arc: results.
20
25
30
Chapter 3
Examples
98
3.1
Examples
Example 1
A frame is clamped clamped at node 1 and has a hinge hinge at node 8. An intern internal al hinge is presen presentt at nodes 2 and 3. At node 4 a point load load of −20 kN is applied and at node 7 a bending b ending momen momentt of −1 kNm is applied. applied. At node 8 a vertical displacement of −0.02m is imposed. The frame consists consists of rectangular rectangular concrete beams (Young’s (Young’s modulus E modulus E = = 30 × 106 kN kN//m2 and Poisson coefficient ν coefficient ν = = 0.2) with a height of 0. 0 .4m and a width of 0. 0 .2 m. 2 3
4
2
3
5 4 6
5
7
9
1
6
y 1z
x
8
Figure 3.1: Example 1: nodes, elements, loads and boundary conditions. % StaBIL StaBIL manual % Ex Exam ampl ple e 1 % Un Unit its s: m, kN % Nod Nodes es=[ =[No NodI dID D Nodes= [1 [1 2 3 4 5 6 7 8 9
X 0 0 0 4 6 6 7 6 2
Y 0 4 4 4 4 3 3 0 2
Z] 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
% ref refere erence nce nod node e
% Che Check ck the nod node e coo coordin rdinate ates s as fol follow lows: s: figure plotnodes(Nodes); plotnodes (Nodes); % Ele Elemen ment t typ types es -> {El {EltTy tTypID pID EltName} EltName} Type Types= s= {1 beam }; b=0.20; h=0.40; % Sections=[SecID A Sections= [1 [1 b*h
ky Inf
kz In I nf
Ixx Iyy Izz] 0 0 b*h^3/12];
Example 1
99
% Mat Materi erials als=[M =[MatI atID D E nu] nu]; ; Materials= [1 30e6 0.2]; % Ele Elemen ments ts=[E =[EltI ltID D Elements= [1 2 3 4 5 6
TypID Typ ID 1 1 1 1 1 1
SecID Sec ID 1 1 1 1 1 1
MatID Mat ID 1 1 1 1 1 1
% con concre crete te n1 1 3 4 5 6 6
n2 2 4 5 6 7 8
n3] 9; 9; 9; 9; 9; 9];
% Che Check ck ele elemen ment t def defini initio tions ns as fol follow lows: s: hold( hold ( on ); plotelem(Nodes,Elements,Types); plotelem (Nodes,Elements,Types); title( title ( Nodes and ele element ments s ); % Deg Degree rees s of fre freedo edom m % Ass Assemb emble le a col column umn matrix containi containing ng all DOF DOFs s at whi which ch sti stiffn ffness ess is % pr pres esen ent t in th the e mo mode del: l: DOF=getdof DOF=getdof(Elements,Types); (Elements,Types); % Re Remo move ve al all l DO DOFs Fs eq equal ual to ze zero ro fr from om th the e ve vect ctor or: : % - 2D analys analysis is: : se selec lect t on only ly UX,UY UX,UY,R ,ROT OTZ Z % - cla clamp no n ode 1 % - re remo move ve the horizo horizont ntal al displ displac acem emen ent t at no node de 8 seldof=[0.03 seldof=[0.03; ; 0.04; 0.05; 1.00; 8.01]; DOF=removedof DOF=removedof(DOF,seldof); (DOF,seldof); % Ass Assemb embly ly of sti stiffn ffness ess matrix matrix K K=asmkm K=asmkm(Nodes,Elements,Types (Nodes,Elements,Types,Sections,Materials ,Sections,Materials,DOF); ,DOF); % No Noda dal l lo load ads: s: -2 -20 0 kN vert vertic ical ally ly on no node de 4 % -1 kNm on node 7. seldof=[4.02 seldof=[4.02; ; 7.06]; PLoad=[-20; PLoad=[-20; -1]; % Ass Assemb embly ly of the load vec vector tors: s: P=nodalvalues P=nodalvalues(DOF,seldof,PLoad); (DOF,seldof,PLoad); % Constr Constraint aint equati equations ons: : Consta Constant= nt=Coef1 Coef1*DOF1 *DOF1+Coef2 +Coef2*DOF2 *DOF2+ + ... % Constra Constraint ints s=[C =[Cons onstan tant t Coe Coef1 f1 DOF1 DOF1 Coe Coef2 f2 DOF2 DOF2 ...] ...] % Imp Impose osed d dis displa placem cement ent -0.02 m ver vertic ticall ally y in nod node e 8. Constr= [0 1 2.01 -1 3.01; 0 1 2.02 -1 3.02; -0.02 1 8 . 02 Na N aN NaN]; % Add con constr strain aint t equ equatio ations ns [K,P]=addconstr [K,P]=addconstr(Constr,DOF,K,P); (Constr,DOF,K,P); % Solve K * U = P U=K\P; % Plot displacements displacements figure plotdisp(Nodes,Elements,Types,DOF,U) plotdisp (Nodes,Elements,Types,DOF,U)
100
% The dis displa placem cement ents s can be dis displa played yed as fol follow lows: s: printdisp(Nodes,DOF,U); printdisp (Nodes,DOF,U); % Com Comput pute e ele elemen ment t for forces ces Forces=elemforces Forces=elemforces(Nodes,Elements,Types (Nodes,Elements,Types,Sections,Materials ,Sections,Materials,DOF,U); ,DOF,U); % Th The e el elem emen ent t fo forc rces es ca can n be di disp spla laye yed d in a or orde derl rly y ta tabl ble: e: printforc(Elements,Forces); printforc (Elements,Forces); % Plo Plot t ele elemen ment t for forces ces figure plotforc( plotforc ( norm ,Nodes,Elements,Types,Forces) title( title ( Normal forces ) figure plotforc( plotforc ( sheary ,Nodes,Elements,Types,Forces) title( title ( Shear forces ) figure plotforc( plotforc ( momz ,Nodes,Elements,Types,Forces) title( title ( Bending moment moments s )
Examples
Example 1
101
−8.823 3.485
−11.177 3.485
3.485 0.000 −11.177 0.000
−8.823
(a) Displacements
(b) Normal forces −11.177
3.485 8.823
3.485
(c) Shear forces
−11.177
−11.177
−3.485 8.823
0.000 0.000
−0.000 −3.485 −3.485 −0.000
−3.485
12.939 12.939 35.293 35.293
−13.939
(d) Bending moments
Figure 3.2: Results of example 1.
−1.000 −1.000 10.454 9.454
0.000
102
3.2
Examples
Example 2 7
y 1z
x
2
1
2
3
3
4
4
5
5
6
Figure 3.3: Example 2: nodes, elements, loads and boundary conditions. % StaBIL manual % Example 2 % Units: m, kN % Nodes=[NodID Nodes= [1 2 3 4 5 6 7
X Y Z] 0 0 0.5 0 1 0 0 -0.25 1 -0.25 1.25 0 0.3 0.4
0; 0; 0; 0; 0; 0; 0];
% reference node
% Check the node coordinates as follows: figure plotnodes(Nodes); % Element types -> {EltTypID EltName} Types= {1 beam ; 2 truss }; b=0.02; h=0.05; A=0.0002; % Sections=[SecID A Sections= [1 b*h 2 A
ky Inf NaN
kz Inf NaN
% Materials=[MatID E nu]; Materials= [1 210e6 0.3];
Ixx Iyy Izz] 0 0 b*h^3/12; NaN NaN NaN];
% steel
% Elements=[EltID TypID SecID MatID n1 n2 n3] Elements= [1 1 1 1 1 2 7; 2 1 1 1 2 3 7;
Example 2
103
3 4 5
2 2 2
2 2 2
1 1 1
1 3 3
4 5 6
NaN; NaN; NaN];
% Check element definitions as follows: hold( on ); plotelem(Nodes,Elements,Types); title( Nodes and elements ); % Degrees of freedom % Assemble a column matrix containing all DOFs at which stiffness is % present in the model: DOF=getdof(Elements,Types); % Remove all DOFs equal to zero from the vector: % - 2D analysis: select only UX,UY,ROTZ % - clamp nodes 4, 5, 6 seldof=[0.03; 0.04; 0.05; 4.00; 5.00; 6.00]; DOF=removedof(DOF,seldof); % Assembly of stiffness matrix K K=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % Distributed loads are specified in the global coordinate system % DLoads=[EltID n1globalX n1globalY n1globalZ ...] DLoads= [1 0 0 0 0 -0.02 0; 2 0 -0.02 0 0 0 0;]; P=elemloads(DLoads,Nodes,Elements,Types,DOF); % Solve K * U = P U=K\P; % Plot displacements figure plotdisp(Nodes,Elements,Types,DOF,U,DLoads,Sections,Materials) title( Displacements ) % The displacements can be displayed as follows: printdisp(Nodes,DOF,U); % Compute element forces Forces=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads); % The element forces can be displayed in an orderly table: printforc(Elements,Forces); % Plot element forces figure plotforc( norm ,Nodes,Elements,Types,Forces,DLoads) title( Normal forces ) figure plotforc( sheary ,Nodes,Elements,Types,Forces,DLoads) title( Shear forces )
104
Examples
figure plotforc( momz ,Nodes,Elements,Types,Forces,DLoads) title( Bending moments )
−5000e−6
0e−6
0e−6
−5000e−6
(a) Displacements
−5000e−6
0e−6
−5000e−6
(b) Normal forces 5000e−6
0e−6
0e−6
−5000e−6
(c) Shear forces
−1667e−6 −1667e−6
(d) Bending moments
Figure 3.4: Results of example 2.
0e−6
0e−6
Example 3
3.3
105
Example 3
The frame consists of rectangular concrete beams (Young’s modulus E = 30 × 106 kN/m2 and Poisson coefficient ν = 0.2) with a height of 0.2m and a width of 0.1 m.
6 8
11 10
12 3 5
15
7
14
7
9 13
2
4 8
6
10
2
3
z y 9 x
5
1
4 1 Figure 3.5: Example 3: nodes, elements, loads and boundary conditions. % StaBIL manual % Example 3 % Units: m, kN % Nodes=[NodID Nodes= [1 2 3 4 5 6 7 8 9 10
X 2.5 2.5 2.5 2.5 1.25 1.25 0 0 0 0
Y 0 5 5 0 0 5 0 5 0 5
Z] 0; 0; 2.5; 2.5; 3.75; 3.75; 2.5; 2.5; 0; 0];
% reference node
% Check the node coordinates as follows: figure
106
Examples
plotnodes(Nodes); % Element types -> {EltTypID EltName} Types= {1 beam }; b=0.1; h=0.2; % Sections=[SecID A Sections= [1 b*h
ky Inf
kz Inf
Ixx 0.229*h/b^3
% Materials=[MatID E nu]; Materials= [1 30e6 0.2]; % Elements=[EltID Elements= [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
TypID 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
SecID 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Iyy h*b^3/12
Izz] b*h^3/12];
% concrete MatID 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
n1 1 2 10 9 1 2 10 9 3 3 8 5 4 7 7
n2 2 10 9 1 4 3 8 7 4 6 6 6 5 5 8
n3] 3; 3; 8; 5; 3; 4; 7; 8; 2; 10; 10; 8; 9; 9; 10];
% Check node and element definitions as follows: hold( on ); plotelem(Nodes,Elements,Types); title( Nodes and elements ); DOF=getdof(Elements,Types); % Remove all DOFs equal to zero from the vector: % clamp nodes 1, 2, 9 and 10 seldof=[1.00; 2.00; 9.00; 10.00]; DOF=removedof(DOF,seldof); % Assembly of stiffness matrix K K=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % Nodal loads: load case 1: 6 kN horizontally on nodes 4, 7, 5 % load case 2: 15 kN horizontally on nodes 7 seldof=[4.02; 7.02; 5.02]; PLoad=[6 0; 6 15; 6 0]; % Assembly of the load vectors: P=nodalvalues(DOF,seldof,PLoad);
Example 3
% Solve K * U = P U=K\P; % Plot displacements figure plotdisp(Nodes,Elements,Types,DOF,U) title( Displacements ) % The displacements can be displayed as follows: printdisp(Nodes,DOF,U(:,1)); printdisp(Nodes,DOF,U(:,2)); % Compute element forces Forces=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U); % The element forces can be displayed in a orderly table: printforc(Elements,Forces(:,:,1)); printforc(Elements,Forces(:,:,2)); % Plot element forces figure plotforc( norm ,Nodes,Elements,Types,Forces(:,:,1)) title( Normal forces ) figure plotforc( sheary ,Nodes,Elements,Types,Forces(:,:,1)) title( Shear forces along the y-axis ) figure plotforc( shearz ,Nodes,Elements,Types,Forces(:,:,1)) title( Shear forces along the z-axis ) figure plotforc( momx ,Nodes,Elements,Types,Forces(:,:,1)) title( Torsional moments ) figure plotforc( momy ,Nodes,Elements,Types,Forces(:,:,1)) title( Bending moments around the y-axis ) figure plotforc( momz ,Nodes,Elements,Types,Forces(:,:,1)) title( Bending moments around the z-axis ) figure plotforc( norm ,Nodes,Elements,Types,Forces(:,:,2)) title( Normal forces ) figure plotforc( sheary ,Nodes,Elements,Types,Forces(:,:,2)) title( Shear forces along the y-axis ) figure plotforc( shearz ,Nodes,Elements,Types,Forces(:,:,2)) title( Shear forces along the z-axis )
107
108
Examples
figure plotforc( momx ,Nodes,Elements,Types,Forces(:,:,2)) title( Torsional moments ) figure plotforc( momy ,Nodes,Elements,Types,Forces(:,:,2)) title( Bending moments around the y-axis ) figure plotforc( momz ,Nodes,Elements,Types,Forces(:,:,2)) title( Bending moments around the z-axis )
Figure 3.6: Displacements for load cases 1 and 2.
Example 3
109
(a) Normal forces
(b) Torsional moments
(c) Shear forces along z
(d) Shear forces along y
(e) Bending moments around y
(f) Bending moments around z
Figure 3.7: Results for load case 1 of example 3.
110
Examples
(a) Normal forces
(b) Torsional moments
(c) Shear forces along z
(d) Shear forces along y
(e) Bending moments around y
(f) Bending moments around z
Figure 3.8: Results for load case 2 of example 3.
Example 4
3.4
111
Example 4
The frame consists of rectangular concrete beams (Young’s modulus E = 35 × 109 N/m2 and Poisson coefficient ν = 0.2) with a height of 0.5m and a width of 0.2 m and concrete columns of 0.3m by 0.2 m.
205 206
307 105
106 203
303 5
3 3
206
203
103
308 106
204
205 305 304 204 6 201 104 210 306 105310 301 309 201 104 202 101 1104
1 10 z y 1 x
6 103
5 101
302 4
2
102
202
102
2 Figure 3.9: Example 4: nodes, elements, loads and boundary conditions. % StaBIL manual % Example 4 % Units: m, kN % Types={EltTypID EltName} Types= {1 beam ; 2 truss }; % Sections bCol=0.2; % Column section width hCol=0.3; % Column section height bBeam=0.2; % Beam section width hBeam=0.5; % Beam section height Atruss=0.001;
112
% Sections=[SecID Sections= [1 2 3
Examples
A hCol*bCol hBeam*bBeam Atruss
% Materials=[MatID E Materials= [1 35e9 2 210e9
ky Inf Inf NaN
nu 0.2 0.3
kz Inf Inf NaN
Ixx 0.196*bCol^3*hCol 0.249*bBeam^3*hBeam NaN
Iyy bCol^3*hCol/12 bBeam^3*hBeam/12 NaN
Izz hCol^3*bCol/12 hBeam^3*bBeam/12 NaN
rho]; 2500; %concrete 7850]; %steel
L=5; H=3.5; B=4; % Nodes=[NodID X Y Z] Nodes= [1 0 0 0; 2 L 0 0] Nodes=reprow(Nodes,1:2,2,[2 0 0 H]) Nodes=[Nodes; 10 2 0 2] Nodes=reprow(Nodes,1:7,2,[100 0 B 0])
% reference node
figure plotnodes(Nodes); % Elements=[EltID TypID SecID MatID Elements=[ 1 1 1 1 2 1 1 1 Elements=reprow(Elements,1:2,1,[2 0 Elements=[ Elements; 5 1 2 1 6 1 2 1 Elements=reprow(Elements,1:6,2,[100 Elements=[ Elements; 301 2 3 2 302 2 3 2 Elements=reprow(Elements,19:20,1,[2 Elements=reprow(Elements,19:22,1,[4 Elements=[ Elements; 309 2 3 2 310 2 3 2
n1 n2 1 3 2 4 0 0 2
n3] 10; 10]; 2 0])
3 4 10; 5 6 10]; 0 0 0 100 100 100]) 3 103 NaN; 4 104 NaN]; 0 0 0 2 2 0]) 0 0 0 100 100 0]) 4 106 NaN; 6 104 NaN];
hold( on ); plotelem(Nodes,Elements,Types); title( Nodes and elements ); % Plot elements in different colors in order to check the section definitions figure plotelem(Nodes,Elements(find(Elements(:,3)==1),:),Types, r ); hold( on ); plotelem(Nodes,Elements(find(Elements(:,3)==2),:),Types, g ); plotelem(Nodes,Elements(find(Elements(:,3)==3),:),Types, b ); title( Elements: sections ) % Degrees of freedom DOF=getdof(Elements,Types);
yt hCol/2 hBeam/2 NaN
y h h N
Example 4
% Boundary conditions: hinges seldof=[ 1.01; 1.02; 1.03; 2.01; 2.02; 2.03; 101.01; 101.02; 101.03; 102.01; 102.02; 102.03; 201.01; 201.02; 201.03; 202.01; 202.02; 202.03;]; DOF=removedof(DOF,seldof); % Assembly of stiffness matrix K K=asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % Loads % Own weight DLoadsOwn=accel([0 0 9.81],Elements,Types,Sections,Materials); % Wind load % DLoads=[EltID n1globalX n1globalY n1globalZ ...] DLoadsWind =[1 0 0 0 0 1500 0; 2 0 0 0 0 1500 0; 3 0 1500 0 0 1500 0; 4 0 1500 0 0 1500 0]; DLoads= multdloads(DLoadsOwn,DLoadsWind); P=elemloads(DLoads,Nodes,Elements,Types,DOF); % Solve K * U = P U=K\P; figure plotdisp(Nodes,Elements,Types,DOF,U(:,1),DLoads(:,:,1),Sections,Materials) title( Displacements: own weight ) figure plotdisp(Nodes,Elements,Types,DOF,U(:,2),DLoads(:,:,2),Sections,Materials) title( Displacements: wind ) % Compute forces [ForcesLCS,ForcesGCS]=elemforces(Nodes,Elements,Types,Sections,Materials,DOF,U,DLoads); % Compute reaction forces for load case 1 Freac=reaction(Elements,ForcesGCS(:,:,1),[1.03; 2.03; 101.03; 102.03; 201.03; 202.03]) % Plot element forces for load case 1 figure plotforc( norm ,Nodes,Elements,Types,ForcesLCS(:,:,1),DLoads(:,:,1)) title( Normal forces: Own weight ) figure plotforc( sheary ,Nodes,Elements,Types,ForcesLCS(:,:,1),DLoads(:,:,1)) title( Shear forces along y: Own weight ) figure plotforc( shearz ,Nodes,Elements,Types,ForcesLCS(:,:,1),DLoads(:,:,1)) title( Shear forces along z: Own weight ) figure plotforc( momx ,Nodes,Elements,Types,ForcesLCS(:,:,1),DLoads(:,:,1))
113
114
title( Torsional moments: Own weight ) figure plotforc( momy ,Nodes,Elements,Types,ForcesLCS(:,:,1),DLoads(:,:,1)) title( Bending moments around y: Own weight ) figure plotforc( momz ,Nodes,Elements,Types,ForcesLCS(:,:,1),DLoads(:,:,1)) title( Bending moments around z: Own weight ) % Plot element forces for load case 2 figure plotforc( norm ,Nodes,Elements,Types,ForcesLCS(:,:,2),DLoads(:,:,2)) title( Normal forces: Wind ) figure plotforc( sheary ,Nodes,Elements,Types,ForcesLCS(:,:,2),DLoads(:,:,2)) title( Shear forces along y: Wind ) figure plotforc( shearz ,Nodes,Elements,Types,ForcesLCS(:,:,2),DLoads(:,:,2)) title( Shear forces along z: Wind ) figure plotforc( momx ,Nodes,Elements,Types,ForcesLCS(:,:,2),DLoads(:,:,2)) title( Torsional moments: Wind ) figure plotforc( momy ,Nodes,Elements,Types,ForcesLCS(:,:,2),DLoads(:,:,2)) title( Bending moments around y : Wind ) figure plotforc( momz ,Nodes,Elements,Types,ForcesLCS(:,:,2),DLoads(:,:,2)) title( Bending moments around z: Wind ) % Load combinations % Safety factors gamma_own=1.35; gamma_wind=1.5; % Combination factors psi_wind=1; % Load combination U_UGT=gamma_own*U(:,1)+gamma_wind*psi_wind*U(:,2); Forces_UGT=gamma_own*ForcesLCS(:,:,1)+gamma_wind*psi_wind*ForcesLCS(:,:,2); DLoads_UGT(:,1)=DLoads(:,1,1) DLoads_UGT(:,2:7)=gamma_own*DLoads(:,2:7,1)+gamma_wind*psi_wind*DLoads(:,2:7,2); figure plotdisp(Nodes,Elements,Types,DOF,U_UGT,DLoads_UGT,Sections,Materials) printdisp(Nodes,DOF,U_UGT); printforc(Elements,Forces_UGT); % Plot element forces figure plotforc( norm ,Nodes,Elements,Types,Forces_UGT,DLoads_UGT) title( Normal forces: UGT ) figure plotforc( sheary ,Nodes,Elements,Types,Forces_UGT,DLoads_UGT)
Examples
Example 4
title( Shear forces along y: UGT ) figure plotforc( shearz ,Nodes,Elements,Types,Forces_UGT,DLoads_UGT) title( Shear forces along z: UGT ) figure plotforc( momx ,Nodes,Elements,Types,Forces_UGT,DLoads_UGT) title( Torsional moments: UGT ) figure plotforc( momy ,Nodes,Elements,Types,Forces_UGT,DLoads_UGT) title( Bending moments around y: UGT ) figure plotforc( momz ,Nodes,Elements,Types,Forces_UGT,DLoads_UGT) title( Bending moments around z: UGT ) % Plot stresses figure plotstress( snorm ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Normal stresses due to normal forces ) figure plotstress( smomyt ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Normal stresses due to bending moments around y: top ) figure plotstress( smomyb ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Normal stresses due to bending moments around y: bottom ) figure plotstress( smomzt ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Normal stresses due to bending moments around z: top ) figure plotstress( smomzb ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Normal stresses due to bending moments around z: bottom ) figure plotstress( smax ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Maximal normal stresses ) figure plotstress( smin ,Nodes,Elements,Types,Sections,Forces_UGT,DLoads_UGT) title( Minimal normal stresses )
115
116
Examples
(a) Displacements
(b) Normal forces
(c) Torsional moments
(d) Shear forces along z
(e) Shear forces along y
(f) Bending moments around y
(g) Bending moments around z
Figure 3.10: Results for load case 1 of example 4.
Example 4
117
(a) Displacements
(b) Normal forces
(c) Torsional moments
(d) Shear forces along z
(e) Shear forces along y
(f) Bending moments around y
(g) Bending moments around z
Figure 3.11: Results for load case 2 of example 4.
118
Examples
(a) Displacements
(b) Normal forces
(c) Torsional moments
(d) Shear forces along z
(e) Shear forces along y
(f) Bending moments around y
(g) Bending moments around z
Figure 3.12: Results for load combination UGT of example 4.
Example 4
119
(a) normal forces
(b) bending moments around y: top
(c) bending moments around y: bottom
(d) bending moments around z: top
(e) bending moments around z: bottom
(f) maximal
(g) minimal
Figure 3.13: Normal stresses for load combination UGT of example 4.
120
3.5
Examples
Example 5
In this example the eigenfrequencies of a simply supported rectangular plate are calculated using shell4 and compared with the theoretical solution. % dynamic plate problem (dkt element) % parameters Lx = 10; Ly = 10; t = 1; E = 200000; nu = 0.3; rho = 7000; m = 20; n = 20; Types = {1 shell4 }; Sections = [1 t]; Materials = [1 E nu rho];
% % % % % % % % % % %
length x-direction length y-direction thickness plate E-modulus Poisson-coeff. mass density number of elements in x-direction number of elements in y-direction {EltTypID EltName} [SecID t] [MatID E nu]
% mesh Line1 Line2 Line3 Line4
= = = =
[0 0 0;Lx 0 0]; [Lx 0 0;Lx Ly 0]; [Lx Ly 0;0 Ly 0]; [0 Ly 0;0 0 0];
[Nodes,Elements,Edge1,Edge2,Edge3,Edge4] = makemesh(Line1,Line2,Line3,Line4,n,m,Types(1,:),1,1); figure; plotnodes(Nodes); figure; plotelem(Nodes,Elements,Types); % DOFs (simply supported plate) DOF = getdof(Elements,Types); sdof = [0.01;0.02;0.06;[Edge1;Edge2;Edge3;Edge4]+0.03;[Edge1;Edge3]+0.05;[Edge2;Edge4]+0.04]; DOF = removedof(DOF,sdof); % K & M [K,M] = asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % eigenmodes nMode = 10; [phi,omega] = eigfem(K,M,nMode); figure; animdisp(Nodes,Elements,Types,DOF,phi(:,1)); % analytical solution [mm,nn] = meshgrid((1:m),(1:n)); aomega = sqrt(E*t^2/(12*(1-nu^2)*rho))*((mm*pi/Lx).^2+(nn*pi/Ly).^2);
Example 5
121
aomega = reshape(aomega,numel(aomega),1); aomega = sort(aomega); ratio = omega./aomega(1:nMode)
(a) 0.1268 Hz
(b) 0.3170 Hz
(c) 0.3170 Hz
(d) 0.5050 Hz
Figure 3.14: The first four eigenmodes of a thin plate. f FEM[Hz] 0.1268 0.3170 0.3170 0.5050 0.6355 0.6355 0.8202 0.8202 1.0848 1.0848
f analytical[Hz] 0.1270 0.3176 0.3176 0.5082 0.6352 0.6352 0.8258 0.8258 1.0799 1.0799
f FEM/f analytical 0.9983 0.9982 0.9982 0.9939 1.0004 1.0004 0.9932 0.9932 1.0046 1.0046
Table 3.1: Eigenfrequencies of a thin plate.
122
Examples
3.6
Example 6
A thin disk with a circular hole is subjected to uniaxial tension. The stress concentration around the hole is examined. % Stress concentration around a circular hole in a disk. r = 1; % Radius of hole [m] L = 20; % Width of plate [m] p = 1; % Uniform pressure [N/m^2] n = 40; % Mesh parameter (nElem = m*n) m = 20; % Mesh parameter Types = {1 shell8 }; % {EltTypID EltName} Sections = [1 1]; % [SecID t] Materials = [1 200000 0]; % [MatID E nu]
% define lines Line1 = [r 0 0;L/2 0 0]; Line2 = [L/2 0 0;L/2 L/2 0;0 L/2 0]; Line3 = [0 L/2 0;0 r 0]; Line4 = [r*sin((0:5). *pi/10) r*cos((0:5). *pi/10) zeros(6,1)]; [Nodes,Elements,Edge1,Edge2,Edge3,Edge4] = makemesh(Line1,Line2,Line3,Line4,... m,n,Types(1,:),Sections(1,1),Materials(1,1), L2method , linear ); % Check mesh: figure; plotnodes(Nodes, numbering , off ); hold( on ) plotelem(Nodes,Elements,Types, numbering , off ); title( Nodes and elements ); % Select all dof: DOF = getdof(Elements,Types); % Apply boundary conditions: % - Line1 and Line3: symmetry condition seldof = [0.03;0.04;0.05;0.06;Edge1+0.02;Edge3+0.01]; DOF = removedof(DOF,seldof); % Assemble K: K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF); % C P U
Apply load to upper edge (equivalent nodal forces): = selectdof(DOF,Edge2(1:(end-1)/2+1)+0.02); = C . *[1/6; repmat([2/3; 1/3],((length(Edge2)-1)/2-2)/2,1);2/3; 1/6]*p*L/m; = K\P;
% Calculate stress in cylindrical global coordinate system [SeGCS]=elemstress(Nodes,Elements,Types,Sections,Materials,DOF,U, gcs , cyl ); % Get nodal stress from element solution [SnGCS] = nodalstress(Nodes,Elements,Types,SeGCS); % plot results figure;
Example 6
123
plotstresscontourf( sx ,Nodes,Elements,Types,SnGCS) title( \sigma_{r} ) figure; plotstresscontourf( sy ,Nodes,Elements,Types,SnGCS) title( \sigma_{\theta} ) figure; plotstresscontourf( sxy ,Nodes,Elements,Types,SnGCS) title( \sigma_{r\theta} )
(a)
σθ
(c)
σrθ
(b)
σr
Figure 3.15: Results.
124
3.7
Examples
Example 7
The geometrical nonlinear response of a truss is examined. clear; %% Geometric nonlinear truss with multiple snap-throughs % based on example on p281 of Non-linear Modeling and Analysis of Solids % and Structures, Krenk (2008) h = 1; % Height truss ll = 1.697; % Lower length hl = 1.414; % Upper length w = 2; % Width Sections = [1 1]; Materials = [1 1 1 1]; Types = {1 truss }; Nodes = [(1:3). ll*(-1:1). zeros(3,2)]; Nodes = reprow(Nodes,(1:3),1,[3 0 w 0]); Nodes = [Nodes; (7:9). hl*(-1:1). (w/2)*ones(3,1) h*ones(3,1)]; Elements = [1 2;2 3;4 5;5 6;1 4;2 5;3 6; 1 7;2 8;3 9;4 7;5 8;6 9; 1 8;3 8;4 8;6 8;7 8;8 9]; nElem = size(Elements,1); Elements = [(1:nElem). ,ones(nElem,3),Elements]; figure; plotnodes(Nodes); hold on; plotelem(Nodes,Elements,Types); hold off; title( Nodes and Elements ); % BOUNDARY CONDITIONS DOF = getdof(Elements,Types); seldof = [(1:6)+0.01,(1:6)+0.02,(1:6)+0.03]. ; DOF = removedof(DOF,seldof); % POINT LOAD seldof = (7:9)+0.03; L = selectdof(DOF,seldof); P = L . *[-1.5;-1;-1.5]; % SELECT VERTICAL DISPLACEMENT MID SPAN seldof = 8.03; Lw = selectdof(DOF,seldof); seldof = 9.03; Lv = selectdof(DOF,seldof); % LINEAR ANALYSIS K = asmkm(Nodes,Elements,Types,Sections,Materials,DOF); U = K\P; % LOAD CONTROL (NEWTON-RAPHSON) Options.method = nr ;
% nonlinear solver: Newton-Raphson
Example 7
125
nloadincrem = ceil(1/0.03); % number of load increments Options.nloadincrem = nloadincrem; [Unr,K,lambda_nr] = nlsolver(Nodes,Elements,Types,Sections,Materials,DOF,P,[],Options); % ARC-LENGTH CONTROL Options.arclength = on ; Options.maxoiter = 3000; % maximum number of arc-length steps Options.initdl = 0.05; % Initial arc-length [Ual,K,lambda_al] = nlsolver(Nodes,Elements,Types,Sections,Materials,DOF,P,[],Options); % PLOT RESULTS figure plotdisp(Nodes,Elements,Types,DOF,U, dispscal ,1); title( Linear small displacements ) figure plotdisp(Nodes,Elements,Types,DOF,Unr(:,end), dispscal ,1); title( Nonlinear large displacements ) v=axis; figure subplot(2,2,1) plotdisp(Nodes,Elements,Types,DOF,[Ual(:,end),Unr(:,end)], dispscal ,1, dispmax , off ); subplot(2,2,2) plot(-Lw*Ual,lambda_al,-Lw*Unr,lambda_nr); axis([-0.5 2.5 -0.1 0.1]);xlabel( w );ylabel( P ); subplot(2,2,3) plot(-Lv*Ual,lambda_al,-Lv*Unr,lambda_nr); axis([-0.5 2.5 -0.1 0.1]);xlabel( v );ylabel( P ); subplot(2,2,4) plot(-Lw*Ual,-Lv*Ual,-Lw*Unr,-Lv*Unr); axis([-0.5 2.5 0 2]);xlabel( w );ylabel( v ); % ANIMATE LOAD DISPLACEMENT CURVE nstep = size(Ual,2); figure for istep = 1:nstep clf subplot(2,2,1) plotdisp(Nodes,Elements,Types,DOF,Ual(:,istep), dispscal ,1, dispmax , off ); axis(v); subplot(2,2,2) plot(-Lw*Ual(:,1:istep),lambda_al(:,1:istep), .- ); axis([-0.5 2.5 -0.1 0.1]);xlabel( w );ylabel( P ); subplot(2,2,3) plot(-Lv*Ual(:,1:istep),lambda_al(:,1:istep), .- ); axis([-0.5 2.5 -0.1 0.1]);xlabel( v );ylabel( P ); subplot(2,2,4) plot(-Lw*Ual(:,1:istep),-Lv*Ual(:,1:istep), .- ); axis([-0.5 2.5 0 2]);xlabel( w );ylabel( v ); pause(0.1) end
126
Examples
7 11
4
18 16
8 5
8 3
12 19 5
14
9 17
9
1 z y 2 x
1
4 13 6
15
6 10 7
2 3
Maximal displacement: 2.85617
(a) Elements
(b) Nonlinear displacements
0.1 0.08 0.06 0.04 0.02 λ
0 −0.02 −0.04 −0.06 −0.08 −0.1 −0.5
0
0.5
1 w
1.5
2
2.5
(c) Vertical load-displacement curve
Figure 3.16: Results.
Chapter 4
Functions — Alphabetical list
accel.tex accel_beam.tex accel_shell4.tex accel_shell8.tex accel_truss.tex addconstr.tex animdisp.tex asmkm.tex b_shell8.tex cdiff.tex coord_beam.tex coord_shell4.tex coord_shell8.tex coord_truss.tex disp_beam.tex disp_shell4.tex disp_shell8.tex disp_truss.tex dispgcs2lcs_beam.tex dispgcs2lcs_truss.tex dof_beam.tex dof_shell4.tex dof_shell8.tex dof_truss.tex eigfem.tex elemdisp.tex elemforces.tex elemloads.tex elempressure.tex elemshellf.tex elemstress.tex fdiagrgcs_beam.tex fdiagrgcs_truss.tex fdiagrlcs_beam.tex
p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p.
4 24 51 61 33 3 14 2 59 19 26 54 62 35 26 55 62 35 27 36 22 48 58 32 18 5 5 4 39 39 39 28 36 28
128
forces_beam.tex forces_truss.tex forceslcs_beam.tex forceslcs_truss.tex gaussq.tex getdof.tex getmovie.tex grid_shell4.tex grid_shell8.tex ke_beam.tex ke_dkt.tex ke_shell4.tex ke_shell8.tex ke_truss.tex kelcs_beam.tex kelcs_shell4.tex kelcs_truss.tex kenl_kbeam.tex kenl_plane4.tex kenl_truss.tex ldlt.tex lincrpoint.tex loads_beam.tex loads_shell4.tex loads_shell8.tex loads_truss.tex loadslcs_beam.tex loadslcs_shell4.tex makemesh.tex msupf.tex msupt.tex nedloadlcs_beam.tex nelcs_beam.tex newmark.tex nlincrpoint.tex nlincrpoint2.tex nlsolver.tex nodalshellf.tex nodalstress.tex nodalvalues.tex patch_shell4.tex patch_shell8.tex plotdisp.tex plotelem.tex plotforc.tex plotlcs.tex plotnodes.tex plotprincstress.tex plotshellfcontour.tex plotshellfcontourf.tex
Functions — Alphabetical list
p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p. p.
24 34 25 34 40 2 15 56 63 22 50 48 58 32 22 49 32 70 69 69 20 67 23 52 61 33 24 52 41 18 18 26 25 19 67 68 66 42 42 4 56 63 11 10 11 12 10 43 43 44