StaBIL-3.0

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


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


47


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


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


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


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


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


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


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


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


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


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


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


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


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


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


’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


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


[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


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



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


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


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


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


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


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


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.


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


’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).


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


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


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


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


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


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



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)


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)


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


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


Section Material UeGCS Options

Section definition [h] Material definition [E nu rho] Displacements (24 * nTimeSteps) Element options {Option1 Option2 ...}

54


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


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


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


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


Pressure

Distributed loads in corner Nodes

[p1localZ; ...] (4 * 1)

Node F See also


[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


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


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


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)


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


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)


.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


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


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;


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.


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);


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


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


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]);


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


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


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


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


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


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];


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


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


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


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


(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


−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


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


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


Example 4

117

(a) Displacements

(b) Normal forces







118

Examples

(a) Displacements

(b) Normal forces






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


(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

StaBIL-3.0

Recommend Documents