2. Solid Mechanics
2.1. Static Analysis
2.1.1. Bolt Stress Analysis
Note
Location of verification case:
verification/solid_mechanics/static/bolt
This case is not included in the unit tests, because there is no comparison value available.
The linear deformations and the stresses of a bolt under external mechanical loading are computed. The example case FE model is a 3dimensional (solid) model consisting of quadratic isoparametric elements, i.e. HE20 hexahedral, TE10 tetrahedral, and PY15 prismatic elements. The model originates from a study performed at the NLR.
The model had been generated with the Patran (c) mesher and converted to B2000 MDL input format. The model exhibits about 150000 degrees of freedom and a global factored matrix of about 400’000’000 nonzero elements. The figure below displays the external surface mesh:
The analysis with B2000++ is controlled with the MDL commands adir
and case
:
adir
case 1
end
case 1
analysis linear
nbc 1
ebc 0
gradients 1
end
The case
option of adir
specifies which case B2000++ should
process (here: process case 1 and the case 1
description in the
case
command specifies the analysis type and other options related
to that particular type of analysis:The analysis
option of
case
indicates the type of analysis to be performed. linear
will launch the B2000++ linear analysis solver., The nbc
option
activates the natural boundary condition set 1
(here: ‘forces’) to
be applied to the current case. In this context the set pertains to
forces applied., The ebc
option specifies the essential boundary
condition set 0
to be applied to the current case. In this context
the set pertains to ‘zero displacement constraints’. The gradients
1
option requires the gradients (strains, stresses) to be computed
and stored on the database.
Solid model viewing capabilities of baspl++ are demonstrated in the
figures below. The plots can be obtained either with the GUI of
baspl++ or with the scripts found in the directory: deformed.py
will produce the deformation plot and stress.py
the rendering of
the von Mises comparison stress on a cut through the semitransparent
solid.
2.1.2. Cantilever Bracket
Note
Location of verification case:
verification/solid_mechanics/static/cantilever_bracket
One of K. J. Bathes incompressible material tests.
2.1.3. Cook Membrane Problem
Note
Location of verification case:
verification/solid_mechanics/static/cook_membrane
The Cook membrane problem is a classical test case for linear static analysis named after the author, R. D. Cook [cook1974] who first reported it. The structure consist of a trapezoidal surface in the xy plane (see figure below and MDL input file for dimensions and material constants).
The structure is clamped along the edges E4 and it is loaded by a edge traction load along the edges E2 with a total load of 1 in the y direction.
The test is run with Q4.S.MITC.E4, Q8.S.MITC, Q9.S.MITC, and T3.S.MITC shell elements, comparing the global ydisplacement at point P3 to the reported solution of 23.91:
2.1.4. Deployable Ring
Note
Location of verification case:
verification/solid_mechanics/static/deployable_ring
The deployable ring problem was proposed by Goto et al. [goto] and it has been studied, among others, by Rebel [rebel], although with shell elements. It possesses the remarkable property that the ring can considerably be reduced in size by twisting it. The accompanying finite rotations in space make the deployable ring a severe test for the finite rotations capabilities of nonlinear beam and shell elements.
The ring has a radius R = 60 and a rectangular cross section shape with a width of 0.6 a height of 6.0 and a modulus of elasticity of 2.0 10^{5} [no units indicated]. It is clamped in one point and a rotation in the radial direction of the point opposite to the clamped point is applied. The rotation ranges from 0 to 2:math:pi.
The test is performed with B2 and B3 elements, the B3 elements obviously giving more accurate results because they approximate the circle better than the B2 elements. The FE mesh consist of a variable number of B2.S.RS or B3.S.RS elements. The figures below display the geometry and mesh points as well as results for a mesh with 36 B2.S.RS elements.
The analysis results with 36 B2.S.RS elements are displayed in the two figures below. The graph displays the applied moment:math:M_xx due to the applied rotation in the global xdirection, i.e the reaction to the applied rotation, as a function of the applied rotation.
2.1.5. Hyperboloid Shell
Note
Location of verification case:
verification/solid_mechanics/static/hyperboloid_shell
A hyperboloid shell model is generated with 3 and 6 node triangular elements for several mesh densities, testing the convergence. The test compares to results from
DoNyun, KlausJurgen Bathe; A Triangular sixnode shell element; Computers and Structures, Volume 87 Issue 2324, December, 2009, Pages 14511460.
A particularity of this test is the generation of the hyperboloid shell: B2000++ epatch meshes cannot generate hyperboloid surfaces. Therefore, a cylinder is generated and the nodes changed after the initial b2ip++ run, see__init__.py.
2.1.6. Lee’s Frame Problem
Note
Location of verification case:
verification/solid_mechanics/static/lee_frame
Lee’s frame is a simple frame structure subject to large displacements
and rotations. The B2000++ model defined in file frame.mdl
follows
the original geometry and loading, and it is clamped at both ends.
Analysis is performed by means of the B2.S.RS ReissnerSimo beam
element. The test also checks the beam local coordinate system
orientation option (defined by a reference node or by a reference
orientation vector).
Two tests are executed, one with applying a force as shown in the
above picture, and another with a ‘displacement loading’
instead. Solver options (see case
command in input file) for the
force case are
analysis nonlinear
increment_control_type hyperplane
step_size_max 10 # To get a curve with some points...
ebc 1 # zero constraints set
nbc 1 # applied force set
join 0 # Needed because more than one epatch
The structure deforms until it reaches an unstable state, after which a new stable state is reached. Using the continuation processor, the complete static solution can be obtained as displayed in the loaddisplacement plot below.
2.1.7. Node Section Forces
Location of case: $B2VERIFICATION/solid_mechanics/static/node_section_forces
This test demonstrates and checks the extraction of beam and rod element section forces of selected elements attached to a node. Element section forces must be extracted when reaction forces at the border of mesh parts are required, such as hinges.
Tests are made with rods and beams, checking the elementlocal to global force (and moment) transformations). The model is a simple truss clamped at x=0 and loaded by a concentrated force (see figure below).
The model is generated with MDL files and a linear analysis is
performed with B2000++. Thereafter the Simples postprocessing scripts
rods.py
and beams.py loop over a list of tuples (nodeid,
elemid1, elemid2, ...)
where nodeid
is the node identifier to
which elements are attached and the elemid
are the rod or cable
element identifiers attached to that node. For each tuple (nodeid,
elemid1, elemid2, ...)
the script loops over the elements and
extracts section forces and moments (if any) of each element
elemid
at the nodenodeid
,transforms the forces and moments to the global coordinate system,
Sums up the resulting forces (and moments) at the node.
In this particular test case the sum of the element forces must be equal to the reaction forces, because the nodes are the constrained nodes.
The test script rods.py
launches the rod test case and is
listed below:
"""Simples script to extract rod forces of selected elements attached
to a node and to sum them up.
"""
import os
import sys
import numpy
import simples as si
numpy.set_printoptions(precision=1)
if os.system('b2000++ rods.mdl'):
sys.exit(1)
verbose = 0
if len(sys.argv) > 1:
verbose = int(sys.argv[1]) # Set to 1 for result output
m = si.Model('rods.b2m')
Branch = m.get_branch(1)
#Extract rod forces (axial forces)
rforces = si.SamplingPointField(m, 'FORCE_ELEMENT', case=1)
Ftot = numpy.zeros(3) # Total force on structure (for equilibrium
# check only)
# Loop over all lists (nodeid, elemid1, elemid2, ...)
for elist in ((1,1,2), (4,3)):
F = numpy.zeros(3) # Sum of forces of selected elements
for eid in elist[1:]:
e = Branch.get_elem(eid)
forcev = numpy.array([rforces[eid][0,0], 0., 0.]) # Local force vector
f = numpy.dot(forcev, e.elbase) # Transform to global
if verbose == 2:
print(e.elbase)
F += f
if verbose == 1:
print(f'Node {elist[0]}: Section forces sum F={F}')
Ftot += F
if verbose == 1:
print(f'Ftot={Ftot}')
if not numpy.allclose(Ftot, numpy.array((0, 1000, 0))):
print(f'***ERRIOR: No equilibrium, Ftot={Ftot}, required=[0,1000,0]')
sys.exit(1)
2.1.8. Raasch Challenge Test
Note
Location of verification case:
verification/solid_mechanics/static/raasch_challenge
The Raasch challenge test is a linear shell test case where drill rotations, i.e shell inplane rotations, play an essential role. The problem has been presented by Knight [knight1997]. The geometry has the form of a clamped ‘hook’, i.e. a thick curved strip with a small radius followed by an arc with a larger radius. The model was used by G. Rebel in his PhD thesis [rebel] to demonstrate the performance of his shell elements in B2000++ However, in the example the MITC shell elements are used.
The FE model consist of two cylindrical patches, both with a length of
20. The first patch has a radius of 46 and spanning an angle of
150 degrees. The second smaller patch has a radius of 14 and spanning
an angle of 60 degrees in the other direction. The mesh can be
parameterized (see the MDL input file raasch.mdl
). Note that, in
contrast to previous versions of this model, both cylindrical mesh
patches are placed in the default branch 1.
is not specified anymore.
The analysis is performed with 4, 8, and 9 node MITC shell elements.
Note that the B2000++ MITC shell elements are 5/6 dof’s per node
elements. To avoid singularities with inplane rotations, B2000++
automatically defines the correct node types: Nodes with constraints
(EBC conditions) and nodes where elements at a given limit angle meet
will become 6 DOF nodes, while all other nodes are 5 DOF nodes or
nodes of another type. The solution (deformed shape and von Mises
stress) for the Q4.S.MITC.E4
elements and a fine grid, with the
initial shape plotted with outline only, is shown below.
The reference solution by Knight for the tip displacement in z direction is 4.9352. Reference solutions for stresses are not given.
2.1.9. Tensile Strip with Hole
Note
Location of verification case:
verification/solid_mechanics/static/strip_with_circular_hole
Stresses and the stress intensity factor in a aluminum tensile strip with a circular hole are computed with solid elements. The test originates from a mesh generated at CIRA s.p.a. for 2 and 3dimensional solid analysis.
The strip with circular hole is described in section Tensile Strip with Hole of the Simples documentation.
2.1.10. Tower Cable
Note
Location of verification case:
verification/solid_mechanics/static/tower_cable
An initially straight cable is loaded by point masses, the cable’s own weight, and point masses. The static equilibrium and the eigenfrequencies around the static equilibrium are computed. This test originates from an older ADINA test case description [bathe1976] and has been converted to the MKS system of units.
The cable is prestressed with an initial stress, initial strain, or
initial force (see input file cable.mdl
). All 3 variants are
tested, checking the final deformed position and the cable forces. The
dead weights are introduced by means of point masses and referenced in
the natural boundary condition case nbc 1
(file cable.mdl
)
with the terrestrial gravitational acceleration in the proper
direction, i.e. the negative yaxis. Note that the vibration analysis
case (case 2) may not reference the natural boundary condition case
nbc 1
, because it would reintroduce the boundary condition a
second time.
The static loaddisplacement response of a node (node A) and the deformed configuration are shown in the figures below. Values correspond to those of the ADINA test case description.
The vibration analysis with prestress results (the equilibrium
position for the load factor 1.0) are listed in the table below. To
perform a restart with a given equilibrium solution the dof_init
option has to be added to the case defining the free vibration analysis
(here: case 2)
case 2
...
dof_init 'DISP.*.10.0.1'
...
end
or in the command line to B2000++
prompt> b2000++ adir dof_init=DISP.*.10.0.1 adir case=2 DBNAME
The first eigenfrequency \(f_1=0.226 Hz\) computed by B2000+++ is identical, to 3 decimal places, to the solution reported in [bathe1976].
A Simples script :file:run.py
computes eigenfrequencies as a
function of the load factor:
"""Compute eigenfrequencies as a function of the load factor
"""
import os
import numpy
import simples as si
from simples import constants as const
def fixmodecycle(db, case, cycle):
# Renames the MODE datsets: B2000++ save the MODE datasets with
# cycle 0 irrespective of the actual cycle for which K is computed!
sets = sorted(db.keys(f'MODE.1.0.0.{case}.*'))
for s in sets:
newname = f"MODE.1.{cycle}.0.{case}.{s.split('.')[5]}"
db[s].name = newname
# Solve nonlinear problem.
os.system("b2000++ define init=2 adir cases=1 cable.mdl")
# Get computed steps
m = si.Model('cable.b2m')
steps = m.get_steps(case=1)
m.close()
x = []
y = []
for c, loadfactor in enumerate(steps):
cycle = c + 1
print(f"Compute eigenmodes for loadfactor {loadfactor}")
os.system(f"b2000++ adir cases=2 adir dof_init='DISP.*.{cycle}.0.1' cable.b2m")
m = si.Model('cable.b2m')
fixmodecycle(m.get_db(), 2, cycle)
result = m.get_eigenfrequencies(case=2, cycle=cycle)
x.append(loadfactor)
y.append(result.frequency)
m.close()
# Create a graph
g = si.Graph2()
g.axes_label_x = 'Load factor'
g.axes_label_y = 'Eigenfrequency [Hz]'
# Add lines
for mode, color in zip((1,2,3), ("RED1", "GREEN1", "BLUE1")):
v = []
for yy in y:
v.append(yy[mode1])
line = g.add_line(x, v, label=f'Mode {mode}')
line.marker = 'o'
line.color = const.RGBCOLORS[color]
g.add_legend()
g.axes_minmax_x = (0., 1.)
g.save('output.svg')
print('Created plot file output.svg')
2.2. Buckling Analysis
2.2.1. Slit Open Annular Plate
Note
Location of verification case:
verification/solid_mechanics/buckling/annular_plate
This test case was proposed by Basar and Ding [basar1992] and is also reported by Buechter and Ramm [buechter1992] and by Rebel [rebel]. It features an annular plate in the xyplane with a slit, with one side of the slit clamped and the other side subject to a line load p in the zdirection.
A nonlinear analysis is performed with the following case definition:
case 1
analysis nonlinear
ebc 1
nbc 1
max_newton_iterations 50
newton_method conventional
max_step 500
step_size_min 0.001
step_size_init 0.01
end
The figure below displays the zdisplacements at points A, B, and C plotted against the applied load for a 32x4 mesh of Q4.S.MITC.E4 elements and the final deformed configuration at the maximum load factor_6.
2.2.2. Linear Buckling of Box Beam
Note
Location of verification case:
verification/solid_mechanics/buckling/box_beam
The linear buckling load of a slender beam subject to axial compression is computed. The beam is clamped on both ends and its section consists of a thinwalled box symmetric with respect to the cross section y and zaxes. Buckling is triggered either by an imposed axial displacement or a force in the axial direction.
Both a beam mesh and a shell mesh are tested:
Beam mesh 1: Linear B2.S.RS elements with boundary conditions imposed on both beam end nodes and section properties imposed by beam properties.
Beam mesh 2: Linear B2.S.RS elements with boundary conditions imposed on both beam end nodes and section properties imposed by a beam section submodel.
The shell mesh consists of element patches of Q elements forming the box. To impose correct boundary conditions for the shell model, i.e. to prevent the end sections from deforming, RBEtype constraints have to be imposed on either beam end, defining the boundary conditions on the reference nodes 100000 and 100001:
The theoretical solution (Euler formula) is readily obtained from the literature:
where \(E\) is the modulus of elasticity, \(J\) the moment of inertia, and \(L\) the length of the clamped beam. With \(E=70 GPa\), \(L=1.0 m\), box width and height \(W=0.02 m\), wall thickness \(t=0.001m\), the critical load is \(P_{cr}= 12671N\).
2.2.3. Channel Section Beam Problem
Note
Location of verification case:
verification/solid_mechanics/buckling/channel_section_beam
A clamped beam with a Csection (‘channel section beam’) described by Betsch et al. [2] is subject to displacement constraints in the global ydirection. The beam is modeled with Q4, Q8, and Q9 shell elements. Results are displayed below for the Q4, Q8, and Q9 meshes.
2.2.4. Curved Panel
Note
Location of verification case:
verification/solid_mechanics/buckling/curved_panel
The test case computes the postbuckling behavior of a slightly curved panel, comparing results with different shell element types.
The panel has a length L=92mm along the straight edges, a radius R=586mm with an angle \(\alpha\) of 9 degrees along the curved edges and a thickness t=1mm. The mesh grid is 7 by 7 nodes in both directions. Due to symmetry only one quarter of the panel is considered.
The material is linear isotropic, with modulus of elasticity E=7*10^{4}N/mm^{2} and poisson=0.3.
Line loads are applied on edge I in the positive z direction. Edge I has suppressed displacements in radial direction and edge IV is simply supported. Edges II and III have symmetry boundary conditions.
No analytical solutions are available for this problem. The results of
the analysis, i.e the buckling load values, can be consulted in the test
file __init__.py
. The figure below display the convergence behaviour
of this problem.
2.2.5. Cylinder under Torsion Load
Note
Location of verification case:
verification/buckling/cylinder_torsion
The torsion buckling load of a cylinder is computed. The cylinder is
clamped on one end and subject to a torsional load on the other free
end. Dimensions of the cylinder as defined by the epatch
command
of the MDL input file:
length 4000
radius 500
phi1 0 phi2 360
thickness 1
The material is linear isotropic (see input file), and the cylinder is meshed with 36 elements (Q4) or 18 elements (Q8/Q9) in circumferential direction and 20 elements (Q4) or 10 elements (Q8/Q9) in axial direction.
Note: The cylinder mesh generator of epatch
does not automatically
close the cylinder if the angles are 0 and 360 respectively. This has
to be done by the user with the join command:
join 0
automatic
end
Analysis is performed with the meshes described above and compared to a reference solution obtained with a fine Q4 mesh 288 by 144 elements. Results are displayed in the table below: While the Q4 and Q8 meshes deliver results with a 15% error with respect to the reference solution, the Q9 mesh performs very well.
Element type 
λ_{1} 
λ_{2} 
λ_{3} 
λ_{4} 
Q4.MITC.E4 ref. 
0.6004 
0.6004 
0.6277 
0.6277 
Q4.MITC.E4 
0.709 
0.709 
0.727 
0.727 
Q8.S.MITC 
0.709 
0.709 
0.727 
0.727 
Q9.S.MITC 
0.6004 
0.6004 
0.6277 
0.6277 
2.2.6. Hinged Cylindrical Shell
Note
Location of verification case:
verification/solid_mechanics/buckling/cylindrical_shell_hinged
A slightly curved cylindrical panel is loaded with a central concentrated outotplane force. Computations are performed with Q4.S.MITC, Q4.S.MITC.E4, Q8.S.MITC, Q9.S.MITC,
Results from:
Ramm E.: ‘The Riks/Wempner approach  An extension of the displacement control method in nonlinear behavior of cylindrical shells’, Recent Advances in Nonlinear Computational Mechanics, Chapter 3. Pineridge Press Limited: Swansea, U.K. , pp 63  86, 1982
Pedro et al.: , ‘A finite strain quadrilateral shell element based on discrete KirchhoffLove constraints’, International Journal for Numerical Methods in Engineering, Vol 64, 2005, pp 1166  1206
2.2.7. Euler Buckling of Beam
Note
Location of verification case:
verification/solid_mechanics/buckling/euler
The test computes the buckling load of a clamped beam modeled by beam elements and subject to compression and to eccentric load application. This test corresponds to one of the classical Euler buckling problems and it illustrates the capability of the B2000++ environment of dealing with multiple parametric runs with a single input file.
The beam has a length of 1m and a cross section of 0.01m by 0.01m, the
cross section being described by a MDL property
command of type
rectangle
:
property 100
beam_section
mid 1
shape rectangle
d1 0.01 # width
d2 0.01 # height
end
The material is a standard isotropic material (steel) described by a MDL
material
command of type isotropic
:
material 1 isotropic
e 210.0e9
nu 0.3
end
The beam is modeled by a MDL epatch
command with 2 node beams:
epatch 1
geometry line
p1 0 0 0
p2 1.0 0 0
eltype B2.S.RS
beam_offsets 0. (EY) (EZ) 0. (EY) (EZ)
beam_orientation 0 1 0
ne1 (NE)
pid 100
end
The beam offsets are parameterized (beam_offset
option): For each run
the EY
, EZ
, and NE
parameters will have to be passed to
B2000++. Several runs with different parameters, together with the
extraction of the buckling load, will define the diagram below.
The analytic solution (Euler buckling) is readily obtained from the
literature, see run.py
script.
The test loop computes the decrease of the buckling load when an eccentric compression force is applied. Each test is executed with the one of the Python utilities, which allows for passing parameters to the input file:
xy=[]
t = 0.01 # Nominal eccentricity (=beam section width)
print('{:10s} {:10s}'.format('Ecc. ratio', 'Pcr/Pcr0'))
for ez in (0., t/100., t/75., t/50., t/25., t/10., t/7.5, t/5., t/2.5, t):
# Run B2000++ with current parameters passed to B2000++
os.system('rm rf beam.b2m')
if si.util.solve('beam.mdl', params={'EY': 0.0, 'EZ': ez}):
model = si.Model('beam')
result = model.get_buckling_load()
model.close()
pcr = result.bloads[0]
if ez == 0.0:
# Pcrit without eccentricity
pcr0 = pcr
# Append to tuple
xy.append((ez/t, pcr/pcr0))
print('{:10.5f} {:10.5f}'.format(ez/t, pcr/pcr0))
else:
print('Run failed for e/t=', ez/t)
2.2.8. LShaped Frame
Note
Location of verification case:
verification/solid_mechanics/buckling/frame_clamped
An Lshape frame is modeled in the xy plane with shell elements. It is clamped on one end and loaded at the other end at (0,0,0) by a force \(F_y\), see figure below. To trigger buckling a small force in the z direction is applied at (0,0,0).
This problem has been reported, among others, by NourOmid and Rankin [3].
To get beyond the buckling load the “hyperplane” increment control
strategy has been selected, see input file b2test.mdl
. The
diagram and the deformed shape plot below illustrate the nonlinear
behavior:
NourOmid and C.C. Rankin: ‘Finite rotation analysis and consistent linearization using projectors’, Computer Methods in Applied Mechanics and Engineering Vol. 93 pp. 353384, 1991
2.2.9. Hexadome
Note
Location of verification case:
verification/solid_mechanics/buckling/healey_hexadome
This test checks the performance of rod elements in the postbuckling range. It consists of a shallow truss structure referred to as a hexadome. The hexadome was first analyzed by Tim Healey [healey]. The rods have a length L=9m, with the height h=1.5m and the rod cross section area A=1cm^{2}. At each of the intersection nodes (where more than two bars are attached to a node), an external force F=10N is applied. The structure is simply supported. Every bar in the structure is modeled with one rod element.
The nodes of the rod elements have three translational degrees of
freedom. For this reason the rod elements can only be simply
supported, the bending stiffness of a rod element being zero. The dome
only buckles as a result of the instability of the construction, not
as a result of the local buckling of a single bar! Due to the forces,
the rods shorten, and at a certain point the dome becomes unstable and
it buckles. Since rods have only stiffness in axial direction only one
element can be used to model a bar over the length L. The material is
linear isotropic, see input file b2test.mdl
.
x
The results displayed below were obtained with increment control type
hyperplane
. To make sure that the whole path is captured, an
initial load factor increment of 0.01, a maximum step size of 0.2 were
selected, up to a final load factor of 5.0. The load factor increases
and decreases several times during the buckling process, as shown in
the loaddisplacement curve for the midpoint node below.
The test checks the final displacement state for a load fact or of 5.0. Intermediate steps are not checked.
2.2.10. Cylindrical Panel (PSC5)
This composite stiffened panel, originating from a Nastran BDF file, is one of the COCOMAT [cocomat] benchmark problems (see [ha1], [ha2], [ha3]). It has the particularity that the shell elements of the stringer feet require an offset with respect to the neutral axis. This can be modeled with B2000++ in two ways
Using ‘air layers’: The air layers are taken into account the calculation of the total volume of the model, resulting in a difference of volume between the two models.
Defining an eccentricity for specific elements.
The panel mesh and boundary conditions, as generated from the converted Nastran BDF file, are visualized with baspl++:
The linear buckling analysis provides critical loads for a number of selected eigenvalues (here: The first 10 lowest eigenvalues):
python3 print_critical_loads.py
Eigenvalues and buckling loads, case 1
Contributing EBC identifiers: 3
Total reaction force F=(2070.46, 34.5649, 118450, amplitude Famp=118468
Mode Lambda Crit. load
2 0.6268 7.426e+04 *
3 0.7199 8.528e+04 *
4 0.7573 8.972e+04 *
5 0.77 9.122e+04 *
6 0.7964 9.435e+04 *
7 0.802 9.501e+04 *
8 0.8079 9.571e+04 *
9 0.8168 9.677e+04 *
10 0.8257 9.782e+04 *
Lambda (eigenvalue): if > 1 no buckling, < 1 buckling, marked by *
Critical load: F * Lambda (if F != 0)
Reaction force: Sum of reactions due to EBC and NBC
A comparison of the lowest eigenmode with an simulation preformed with Abaqus 6.4:
The static solution up to a load factor of 3, i.e. an end shortening of 3mm, is obtained in a few steps, as can be seen in the load step vs edge shortening graph below:
Again, the panel deformed shape of the last step is obtained with baspl++:
The implicit dynamic (transient) solution, obtained with the Nordsieck solver, gives the dynamic buckling solution  a more accurate solution, exhibiting intermediate buckling states:
2.2.11. Cylindrical Panel (Verolme)
Note
Location of b case:
verification/solid_mechanics/buckling/verolme_panel
This test case is a variant of the panel described by Verolme [verolme1995] and Remmers [remmers1998]. The panel illustrates the postbuckling behaviour of an initially undeformed panel, which ‘jumps’ from its prebuckling state to a stable postbuckling state. Remmers [remmers1998] simulates the behavior by preforming a dynamic (transient) analysis at a load just after the instability point, following the modejumping technique according to Riks and Rankin [riks1996], [riks1997].
The term mode jumping is often used to describe a sudden dynamic change in a static process. In computational mechanics, the abrupt change in wave numbers (such as from pre to postbuckling state) is indicated as a modejump. This phenomenon was first mentioned by Stein [stein1959], describing a buckling experiment with a flat panel. Modejumping is an important phenomenon for stiffened panels, where a modejump may occur from local skin buckling (between the stiffeners) to a global buckling pattern.
In order to simulate such jumps with FEA, one can use a transient solving routine including some form of damping, such as Rayleigh damping, in order to find a stable static path [2]. The artificial damping scheme of B2000++ allows for finding a stable path with ‘static’ analysis. This test case uses this artificial damping method to find a stable postbuckling solution (a nondamped simulation will fail).
The panel is made out of 2024T3 aluminum and simply supported at the straight edges and clamped at the curved edges, as shown in the figure below. The panel is compressed at one of the curved edges, creating an edge shortening of 3 mm.
The geometry and material parameters are listed below:
Panel length \(L=0.5 m\), \(W=0.356 m\)
Panel radius \(R=0.38 m\), panel thickness \(t=0.001 m\)
Young’s modulus \(E=73.1 GPa\), Poisson ratio \(\nu=0.33\)
The numerical model for the test case uses a mesh of 20x20 elements. This will execute quickly, but is too coarse for an accurate result. An 80x80 mesh will produce an accurate result. Higher mesh densities do not show any significant improvement.
The following options are defined in the cases
command:
cases
case 1
residue_function_type artificial_damping
dissipated_energy_fraction 2.e6
step_size_init 0.01
step_size_max 0.005
step_size_min 1e13
max_newton_iterations 30
max_divergences 1
increment_control_type load
end
end
The figure below shows the deformed panel and the loadshortening curve of the analysis for a 160x160 mesh and the loadshortening curves for other mesh densities.
2.3. Transient (Dynamic) Analysis
2.3.1. Forced Vibration (Bar)
Note
Location of verification case:
verification/solid_mechanics/dynamic/forced_vibration_bar
This case tests both the sfunction
and the stabulated
options
for prescribing timedependent natural boundary conditions nbc
,
i.e for the present case, timedependent excitation forces.
A single DOF system is subject to a sinusoidal excitation force in the
global xdirection in the geometrically linear range. The theoretical
solution can be found in the text book by Clough and Penzien and is
displayed in the figure below. The code is contained in the test file
__init__.py
.
The corresponding FE problem is formulated with 1 or several rod
elements in the xdirection in a line. All DOF’s except DOF 1 are
removed for all nodes. The mass is concentrated at the end node and
introduced with a PMASS3 element. The linear dynamic solver is
selected, since the problem is linear. Setting the
multistep_integration_order
to 2
is the sufficient, the
numerical damping being dependent on the time step only. Check with
the input file r2_tabulated.mdl
for the definition of the time
step.
The test checks forced vibration defined by the the case option
nbc
, sfunction
defining the excitation function with an
expression
nbc 1 sfunction 'sin(6.6666666666667*t)'
or stabulated
with pairs (time, value):
nbc 1 stabulated [ 0.0 0.0 ... 2.4740042147 0.707095118523 ]
Both options define the same sinusoidal excitation function of the natural boundary condition. The response is displayed below:
2.3.2. Nonlinear 1 DOF system
Note
Location of verification case:
verification/solid_mechanics/dynamic/one_dof_nonlinear
This simple rod snap problem consists of two rods. The FE model is show below. The rod dimensions are B=H=10 and the sections are circular with a radius of 0.5641896, i.e. a cross section area of 1.0. The material is linear isotropic with a modulus of elasticity E=1.e6. The rods are restrained at nodes 1 and 3, node 2 being free to move in the x and ydirection, and the node 2 is initially deformed to 1 at time 0 and released. The example has been described by Argyris et.al [1], where the analytical formulation of the equation of motion of the two rod system can be found.
The figure below plots the FE solution (node 2, dof 2) vs the analytical one:
Nonlinear oscillations using the finite element technique; J.H. Argyris, P.C. Dunne, T. Angelopoulos; ISD Report 36, Stuttgart, (1972).
2.3.3. Cube in Gravitational Field
Note
Location of verification case:
verification/solid_mechanics/dynamic/one_element_gravity
A solid cube is placed in a gravitational field g=(0, 0, 9.81) at time t=0 and released. The cube must move in ‘free fall’ along the negative zaxis with
and with the falling speed \(v = g_zt\), the acceleration \(a = g_z\) being constant. The model consists of HE8.S, HE20.S.TL, or HE27.S.TL elements. Different mesh densities are tested. The figure below displays the displacement, velocity, and acceleration functions for an integration time of 2s (solid lines are theoretical values).
The dynamic solution is relatively sensitive with respect to the error
tolerance. Relevant tolerances of the case
option (see input file
he.mdl
)
tol_residuum 2e3
tol_solution 1e4
have been found to lead all element types to convergence.
2.4. Free Vibration
2.4.1. Cylindrical Panel
Note
Location of verification case:
verification/solid_mechanics/freevibration/cylindrical_panel_1
The eigenfrequencies and eigenmodes of a thin isotropic cylindrical panel are calculated. The purpose of the analysis is to examine the performance of the MITC shell elements combined to point mass elements (PMASS).
The panel is supported with four springs attached to the corners of the plate. Thus, the panel is not fully constrained and a shift of 1.0 has to be applied in order to be able to solve the eigenvalue problem.
The panel has a length a = 1.039_m, radius R = 1. m, an angle φ = 49.4 and a thickness t = 0. m. The springs have a height h = and a cross section area of 1. m^{2}. The materials are linear isotropic and the properties for the springs and the panel are listed below.
Panel: E=74e9 MPa, Poisson=0.3, density=2700 kg/m^{3}
Springs: E=730 MPa, Poisson=0.3, density=0.024 kg/m^{3} 0. 0.024
The free vibration analysis reveals that the first 6 eigenfrequencies are due to the vibration of the panel on its foundation. Therefore, only the 7^{th} to 10^{th} eigenfrequencies are listed in the table below:
Element 
Mode 7 
Mode 8 
Mode 9 
Q4.S.MITC 
9.94 
12.9 
22.6 
Q8.S.MITC 
9.93 
12.7 
22.4 
Q9.S.MITC 
9.91 
12.7 
22.4 
2.4.2. Eigenfrequencies of Hinged Plate
Note
Location of verification case:
verification/solid_mechanics/freevibration/hinged_plate
The eigenfrequencies and eigenmodes of a thin flat isotropic plate
with a hinge are calculated. The purpose of the analysis is to
demonstrate and check the modelling of a hinge, by means of different
patches and the LINC
command. Mass matrices are calculated with
the default B2000++ consistent mass configuration.
The plate has a total length of a = 1.0 [m], a width b = 0.4 [m] and a thickness t = 0.004 [m]. The material of the plate is linear isotropic with a modulus of elasticity E = 74.E9 [N/m:sup:2], a Poisson number ν = 0.3 and a density ρ = 2700 [kg/m:sup:3]. One of the short edges is clamped, while the others are free. The hinge is located halfway between the two smaller edges (at a = a/2). The input options are listed below.

Selects the hinge type: htype=1 adds a hinge, htype=2 means no hinge. 
The results for the different element types for the plate with or without hinge are listed in the table below:
Element 
Hinge 
f1 
f2 
f3 
f4 
f5 
Q4.S.MI TC.E4 
yes 
3.46 
17.84 
21.56 
56.76 
60.77 
Q4.S.MI TC.E4 
no 
9.73 
17.79 
54.88 
60.75 
84.20 
2.4.3. Eigenfrequencies of Flat Isotropic Plate
Location of verification case: $B2VERIFICATION/solid_mechanics/freevibration/isotropic_plate
The eigenfrequencies and eigenmodes of a thin flat rectangular isotropic plate are calculated for MITC shell elements with 4, 8, and 9 nodes. The plate has a length a = 10 m, width b = 0.4 m and a thickness t = 0.002 m. The material is linear isotropic, with modulus of elasticity E = 72 MPa, Poisson number 0.3, and density 2400 [kg/m^{3}]. The edges I and III are simply supported. The other two edges are free. The plate is modeled using a regular 10 by 10 element mesh.
Analytical solutions can be obtained from various sources, and
computed values are not listed here (see __init__.py
). The
figure below shows the second, third, and fourth eigenmode of the
plate with respect to the initial shape.
f_{1} 
f_{2} 
f_{3} 
f_{4} 
5.0 
17.1 
20.2 
38.6 
2.4.4. Composite Panel with Stringers
Note
Location of verification case:
verification/solid_mechanics/free_vibration/panel3
The three string composite panel is a test case for demonstrating the B2000++ capabilities of dealing with lightweight composite structures modeled with shell or solid elements. A cylindrical panel with 3 stiffeners is subjected to axial compression load. The free vibration modes are to be determined under zero compression load.
The boundary conditions are displayed in the figure below: Both curved
edges are clamped and both straight edges are free. One of the curved
edges is pushed towards the panel, causing it to buckle. All boundary
conditions are essential boundary conditions (ebc
command in the
MDL input files).
The shell model consists of Q9 MITC shell elements. The cylindrical shell and the shell plus the stringer foot are modeled with one element through the thickness, as is the stringer. The laminate orientations and laminate placement are visualized with baspl++:
Free vibration modes and their eigenfrequencies of the panel without prestress of the shell model are computed for load case 1, requiring the first 10 modes:
case 1
ebc 1
analysis free_vibration
nmodes 10
gradients 1
end
The analysis of the shell composite model is launched by the shell command
prompt> b2000++ adir 'cases=1' l "info of solver in cerr" shell.mdl
Eigenfrequencies can be extracted and printed with
b2print_eigenfrequencies and plotted with the baspl++
script view_vibration_modes.py
.
b2print_eigenfrequencies shell
Eigenfrequencies, case 1
Mode FREQ OMEGA MODK MODM
1 205.3 1290 116.1 6.979e05
2 268.7 1688 339.5 0.0001191
3 304.3 1912 497.9 0.0001362
4 308.4 1938 364.3 9.703e05
5 349.6 2197 425.2 8.813e05
6 380.1 2388 296.5 5.198e05
7 434.7 2731 774.1 0.0001038
8 508.5 3195 448.3 4.392e05
9 547.4 3439 414.3 3.503e05
10 606.4 3810 1027 7.073e05
Similarly, the analysis of the solid composite model is launched by the shell command
prompt> b2000++ adir 'cases=1' l "info of solver in cerr" solid.mdl
Eigenfrequencies can be extracted and printed with
b2print_eigenfrequencies and can be plotted with the
baspl++ script view_vibration_modes.py
.
b2print_eigenfrequencies solid
Eigenfrequencies, case 1
Mode Frequency Omega Modal_K Modal_M
1 205.5 1291 7.086e05 118.1
2 265.5 1668 0.0001237 344.4
3 303.5 1907 0.0001339 486.9
4 305.6 1920 8.978e05 331.1
5 351.4 2208 8.546e05 416.6
6 383 2407 5.331e05 308.8
7 439.7 2763 0.000104 793.9
8 513.2 3225 4.565e05 474.7
9 554.8 3486 3.647e05 443.1
10 613.7 3856 7.478e05 1112
2.4.5. PMASS free vibration tests
Note
Location of verification case:
verification/solid_mechanics/freevibration/pmass1
This free vibration tests check the PMASS3 point mass elements. The test compute the first n eigenfrequencies of a clamped beam with a consistent mass matrix or with equivalent PMASS3 elements:
beam.mdl
: Consistent beam element mass formulation.pmass3.mdl
: PMASS3 point masses lumped at the mesh nodes. The density of the beam elements is then set to 0. Note that the PMASS tests are performed with B2 elements only, the B3 elements not being suited for this type of test because of the extra node.
The beam dimensions can be consulted in the params.mdl
file. The
material is isotropic, and the beam is modeled with B2.S.RS elements.
The test compares all configurations, which must give the same result.
2.4.6. PMASS6 free vibration tests
Note
Location of verification case:
verification/solid_mechanics/freevibration/pmass2
This free vibration tests check the PMASS6 point mass elements. The tests compute the first n eigenfrequencies of a clamped beam with various configurations of PMASS6, or equivalent RBE elements. The beam is massless and acts as a spring for the PMASS6 and RBE elements.
Clamped beam with PMASS3, PMASS6, or RBE element at the tip.
Clamped beam with a PMASS6 or RBE element at the tip with an offset y in the xy plane (DOF’s 3 4 5 eliminated) and with an offset in the xz plane (DOF’s 2 4 6 eliminated).
Clamped beam with a PMASS3 element at the tip with an offset in y simulated by a rigid beam.
The beam dimensions can be consulted in the params.mdl
file. The
material is isotropic, and the beam is modeled with B2.S.RS elements.
The mass is a cube of dimensions LX*LY*LZ and density RHO2 (see
params.mdl
) and is modeled with a single HE8 elemen (for the RBE
tests).
The test compares equivalent configurations, i.e. PMASS6 vs. RBE without offset and PMASS6 vs. RBE with and offset.
2.4.7. Prestressed Beams
Note
Location of verification case:
verification/solid_mechanics/freevibration/prestressed_beams
The eigenfrequencies and eigenmodes of a simplified helicopter rotor modeled with beam elements are calculated. The rotor consists of four blades connected to each other by springs.
The geometry data as well as the material properties (all isotropic elastic) are listed below:
\(length L=1.7m\) \(area A=0.01m^2\) \(head width a=0.2m\) \(shaft height h=0.1m\) \(Iyy=Izz=2.8e8m^4\) \(Ip=2*Iyy\) \(Sy=Sz=100m^2\)
\(E_rotor=1.E10{N/m^2}\) \(\rho_rotor=117.6{kg/m^3}\) \(E_springs=1.E6{N/m^2}\) \(\rho_springs=100{kg/m^3}\)
The end point of the shaft is fully clamped.
he solutions are independent of the mesh type or join type and vary only due to the prestress. The table below lists the 10 lowest eigenfrequencies.
Prestress 
1st 
2nd 
3rd 
4th 
5th 
0.0 
2.987 
2.987 
2.987 
2.987 
2.987 
1E7 
44.178 
44.189 
44.189 
44.189 
44.189 
Prestress 
6th 
7th 
8th 
9th 
10th 
0.0 
3.1136 
3.1138 
3.197 
18.643 
18.643 
1E7 
44.3967 
44.3972 
44.590 
130.31 
130.35 
The input options to the test program are:
meshtype: Type of mesh to be generated: meshtype=1 selects one branch, meshtype=2 selects four branches.
jointype: Branch join type: jointype=1 selects automatic join, jointype=2 selects manual join (meshtype=2 only).
prestr: Specifies the value of the prestress: Analysis is performed with prestr=1.0e7) or without prestress (prestr=0).
2.5. Material Failure
2.5.1. Hashin Criterion Failure Envelope
Note
Location of verification case:
verification/solid_mechanics/failure/hashin
The Hashin criterion [hashin80] failure envelope for \(\sigma_{22}\) vs \(\sigma_{12}\) is computed with one single HE solid element by executing B2000++ runs with displacement constraints corresponding to the desired stress state. The results are compared to data from literature [pinho].
Note: The Hashin failure criterion data are entered with the
failure
subcommand of the MDL material
command. Example from
test:
failure hashin
t1 1000. # (not used)
c1 1000. # (not used)
t2 37.5
c2 130.5
s12 66.5
s23 49.09 # equation 68 in [pinho] = 0.378 * c2
end
2.5.2. Solid Elements Rod
Note
Location of example case:
verification/solid_mechanics/kinematic_coupling/pulled_rod_solid_elements
A rod meshed with HE8 solid elements and composed of two sections with two different square cross section areas is clamped at one end and pulled at the other end. The two sections are coupled with the common mesh refinement kinematic coupling method of B2000++. The analysis is linear. Square areas with an area ration of 5:1 and 2.5:1 are checked, the first section area always being the same and the second one being varied.
The analytic reference solution for the displacement \(u_x\) at \(x=2L\) is simple, since the model is ‘statically determined’:
Displacements in the xdirection agree within 0.01% with the analytic solution.
2.6. Model Reduction
2.6.1. Vibration of a Clamped Beam
Note
Location of verification case:
verification/solid_mechanics/model_reduction/beam
The test computes vibration modes of a clamped beam with the Component Mode Synthesis method (also referred to as the CraigBampton method) globallocal analysis method. The analysis consists of 2 steps, one with a global coarse mesh, and one with a local fine mesh. In a first step the global model is reduced to a macroelement. Then, the local model uses this reduced system defined by the global model for making a local (detailed) vibration analysis.
The globallocal analysis is performed as follows:
Condensate the global model
x_global.mdl
.Create the input database of the local model
detail.mdl
.Update the input database of the local model with the information from the global model. This step has to be performed with a Python script (here:
b2ipupdate.py
).Compute vibration modes with the local model.