Examples

The B2000++ examples collection illustrate typical B2000++ use cases. All examples are explained step by step and can be reproduced by the user. Examples cases are an integral part of the verification suite. The examples are distributed with B2000++ and are located under

<prefix>/share/b2000pp

The examples also constitute test cases for the verification of B2000++. And they serve as additional checks for Simples and the MemCom data manager via the b2testrunner unit test verification program.

All analysis and post-processing steps of the cases presented in the examples here are launched by means of scripts. Specific data processing capabilities by the Simples and the baspl++ post-processor are described in the text and also documented in the scripts.

All examples are included in the B2000++ test procedure.

The B2000++ examples must be executed with the Python 3 version of B2000++. Please make sure to activate the correct Python 3 module file.

B2000++ examples can be run one by one (see relevant README files) or all together with the b2testrunner program.

Examples cases are executed with the b2testrunner program:

b2testrunner .

To redirect output to a file :

b2testrunner  . >b2testrunner.out 2>&1

Run n tests concurrently with the -j option. Example: Run with 12 threads

b2testrunner -j 12 .

Solid Mechanics: Buckling

Collapse Analysis of Cylinder with Cutout

Note

Location of example case:

$B2EXAMPLES/buckling/cylinder_cutout

A cylinder with two opposite square cutouts is clamped at one end and compressed at the other end in the axial direction by an imposed displacement of 0.00425. This problem has been proposed by Almroth and Holmes [almroth1972] and has later been taken up by Knight and Rankin [knight2006]. The goal of this example is to demonstrate the collapse analysis and the post-processing capabilities of B2000++. Note that the B2000++ Verification Manual treats this case, too, but with Q4 and Q9 elements.

Similarly to Knight and Rankin, the problem is modeled with an eight of a cylinder (half in the longitudinal direction and a fourth in the circumferential direction):

Cylinder with square cutouts: Mesh. 1/8 of cylinder is modeled.

To this end, 2 cylindrical epatch patches are defined. The figure below displays the mesh with Q4.S.MITC elements. Boundary conditions are symmetry along both longitudinal edges (edges defined by blue points P1-P4 and green points P2-P3) and along one circumferential edge (edges defined by blue points P3-P4), the other circumferential edge being clamped.

EPATCHes and mesh.

Boundary conditions and local node coordinate systems.

The test is run with Q4.S.MITC.E4 shell elements, comparing the cutout tip displacement-load function (point A) against the reported solution by Knight and Rankin.

The analysis control block contains a standard set of options, the non-linear analysis method being the default continuation method):

case 1
   analysis nonlinear
   ebc 1
   ebc 2
   step_size_init 0.05
   step_size_max 0.05
   gradients 1
   rcfo_restrict nodeset "EndShortening"
end

The maximum step size limitation is imposed in order to obtain a load curve with some points. Note the rcfo_restrict option: It contains the node set(s) containing all nodes involved in the imposed displacement, thus rendering post-processing very easy (see stroke-load diagram below).

The y-axis of both diagrams below are scaled to the reaction force \(P/P_{0}\), with \(P_{0} = {{1.2}\pi Et^{2}}\), where \(E=1 \times 10^{7}\) and \(t=0.014\).

The Simples script run.py (1) runs the analysis for selected mesh densities, extracting the graphs for each analysis and (2) collects the graphs and produces the stroke versus load disp_load.png and stroke_load.png out of plane versus load plots.

# Python script for executing a series of B2000++ runs with increasing
# mesh density for different element types, followed by a synthesis
# graph.

import os

import simples as si

# Comparison values obtained with STAGS
stags_dx = [0.00891496, 0.0137947, 0.0204575, 0.0290596, 0.0404457,
            0.0551476, 0.0731026, 0.0914018, 0.0975327, 0.0991906, 0.0989717,
            0.102444, 0.105009, 0.107511, 0.109607, 0.112328]
stags_load = [0.086217, 0.104633, 0.119883, 0.137595, 0.162698,
              0.205865, 0.277067, 0.365396, 0.373607, 0.375367,
              0.379707, 0.386628, 0.394956, 0.395191, 0.38393,
              0.376657]

os.system("rm -rf *.b2m *.modeler")

# Compute P0 = Analytical buckling load of cylinder
E = 1e7
t = 0.014
P0 = 0.25*1.2*3.14159*E*t**2

# Run analysis for different mesh densities for one element type
eltype = 'Q4.S.MITC.E4'
radial_disp_reaction_graphs = []
stroke_reaction_graphs = []
meshes = ((18,32), (36,64))
for ne in meshes:
    mesd = '%s %dx%d' % (eltype, ne[0], ne[1])
    print('Run for mesh density %s' % mesd)
    # Run B2000++ for current parameters
    if si.util.solve('cutout.mdl', params={'ELTYPE': eltype, 'NE1': ne[0], 'NE2': ne[1]}):

        # Open db and get node id of point A.
        model = si.Model('cutout', 'or')

        # Get internal node number of epatch 1 point P3 == point A.
        point = model.get_node_set('EPATCH-1-P')[2]

        # Get all analysis steps values. Scale because TIME goes from
        # 0 to 1, but displacement constrains from 0 to 0.1
        steps = model.get_steps(scale = 0.1)

        # Get the displacements of node *point*. Extract the
        # node-local component 'radial', which is dof 1.
        dy = model.get_dof_value_of_steps('displacements', node = point, dof = 1, system = 'nl')

        # Get sum of all reaction forces due to displacement loading and scale.
        ry = model.get_reaction_force_of_steps(scale = 1./P0)

        model.close()

        # Create Graph for current model and add to list
        radial_disp_reaction_graphs.append((mesd, (dy, ry)))
        stroke_reaction_graphs.append((mesd, (steps, ry)))

# Plot all unctions
colors = ('red', 'green', 'blue', 'cyan', 'magenta', 'yellow')
markerw = 0.5
markers = 3.

g = si.Graph2()
for i, a in enumerate(radial_disp_reaction_graphs):
    xy = a[1]
    line = g.add_line(xy[0], xy[1])
    line.marker= 'o'
    line.markeredgewidth = markerw
    line.markersize = markers
    line.markerfacecolor = None
    line.color = colors[i]
    line.label = a[0]
line = g.add_line(stags_dx, stags_load)
line.color = 'black'
line.label = 'STAGS'
g.axes_label_x = 'Radial displacement at A'
g.axes_label_y = 'Load P/P0'
g.add_legend((.2, .7))
g.save('disp_load.png')

g = si.Graph2()
for i, a in enumerate(stroke_reaction_graphs):
    xy = a[1]
    line = g.add_line(xy[0], xy[1])
    line.marker= 'o'
    line.markeredgewidth = markerw
    line.markersize = markers
    line.markerfacecolor = None
    line.color = colors[i]
    line.label = a[0]
g.axes_label_x = 'Stroke'
g.axes_label_y = 'Load P/P0'
g.add_legend((.2, .7))
g.save('stroke_load.png')

Stroke versus load diagram (Q4 mesh).

Out-of-plane versus load diagram (Q4 mesh).

Buckling of Three Stringer Composite Panel

Note

Location of example case:

$B2EXAMPLES/buckling/panel_3stringers

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 lowest buckling loads are to be determined.

Three Stringer Composite Panel: Dimensions and laminate setup.

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) Thus, the resulting reaction force due to the edge displacement in axial direction is equal to the buckling load when multiplied by the computed eigenvalue.

Three Stringer Composite Panel: Boundary conditions at nodes (red dots).

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++:

Three Stringer Composite Panel (solid model) : Laminates and material axis x (red) and y (green). Laminate 1 is material 11 (blue), laminate 2 is material 12 (yellow), and laminate 3 is material 13 (white).

Buckling modes without prestress of the shell model are computed for load case 1, requiring the first 10 modes:

case 1
  ebc       1
  analysis  linearised_prebuckling
  nmodes    10
  gradients 1
end

The buckling analysis is launched by the shell command

prompt> b2000++ --adir 'cases=1' -l "info of solver in cerr" shell.mdl

Linear buckling analysis results can be extracted and printed with b2print_buckling_loads:

b2print_buckling_loads shell

Linear reaction forces +2.9477108e+05

Number of computed buckling modes for case 1: 10
Mode Load factor    Critical force
   1 +1.7227278e-01 +5.0781032e+04
   2 +1.7237261e-01 +5.0810459e+04
   3 +1.7464288e-01 +5.1479669e+04
   4 +1.7469178e-01 +5.1494083e+04
   5 +1.7908786e-01 +5.2789922e+04
   6 +1.7965825e-01 +5.2958056e+04
   7 +1.8222893e-01 +5.3715819e+04
   8 +1.8233669e-01 +5.3747582e+04
   9 +1.8770931e-01 +5.5331277e+04
  10 +1.8883024e-01 +5.5661694e+04

The first buckling mode shape plot is obtained with the baspl++ script view_buckling_modes.py.

Three Stringer Composite Panel (shell model): First buckling mode.

Buckling of Composite Plate

Note

Location of example case:

$B2EXAMPLES/buckling/plate_buckling

A simply-supported thin quadratic laminated plate of 500 mm by 500 mm is loaded with a forced displacement along one of the plate edges. The linearized buckling load as well as the post-buckling reponse due to the forced displacement of the edge is to be determined.

The plate is simply supported on all four edges, one edge being free to move in the x-direction to allow for the forced displacement. The plate is loaded with a forced displacement loading condition along the free edge 2.

Geometry and loading of plate

The plate is made of a composite material, 16 layers of 0.125mm thickness, with the symmetric stacking sequence (in degrees) \((+45,0,-45,0,+4,0,-45,+90)_s\). The material data are listed below:

\(E_1 = 146.86\ 10^{3} N/mm^2\)

\(E_2 = E_3 = 9.65\ 10\ 10^{3} N/mm^2\)

\(\nu_232 = \nu_13 = 0.023\)

\(G_12 = G_23 = G_13 = 45.5\ 10^{2} N/mm^2\)

FE Model

Now we explain in detail how to get to the solution. We start by defining the Finite Element model with the B2000++ model description language (MDL). MDL contains simple mesh generation capabilities with element patches. A single epatch defines the regular plate mesh and links to the material.

(et="Q4.S.MITC.E4")
(ne=40)

# Mesh definition
epatch 1
  geometry plate  # Geometry type
  p1 0.   0.   0. # 4 plate vertices
  p2 500. 0.   0.
  p3 500. 500. 0.
  p4 0.   500. 0.
  eltype (et)     # Element type
  ne1 (ne)        # N. of elements in first direction
  ne2 (ne)        # N. of elements in second direction
  mid 2           # Link to material identifier
end

Note the parametrization of the element type and the number of elements. Here we select Q4 four node quadrilateral MITC elements.

The plate geometry is uniquely defined by the four corner points p1, p2, p3, and p4 specified counter-clockwise, defining a flat plate. Note that the plate thickness is not needed if the material is a laminate material.

A laminate is defined with the MDL material block, by listing all layers one after the other, each layer specifying the thickness, the layer orientation, and the material identifier.

# Composite stacking sequence definition
material 2
  type laminate
  0.125  45 1
  0.125   0 1
  ...
end

For our model the default global material orientation is selected, i.e the material angle is defined with respect to the global coordinate system as displayed in the figure above.

The composite material constants are defined with a MDL material block, specifying the orthotropic material type.

material 1
  type orthotropic
  e1 146.86e+03
  e2 9.65e+03
  e3 9.65e+03
  nu12 0.3
  nu13 0.3
  nu23 0.023
  g1 45.5e+02
  g2 45.5e+02
  g3 45.5e+02
  failure max_principal_strain
    t +4500.e-6 # tension, micro-strain
    c -3000.e-6 # compression, micro-strain
    s +5000.e-6 # shear, micro-strain
    filter max_of_element
  end
end

A simple failure criterion is used, based on the maximum of the principal eigenvalues of the strain tensor. Its purpose is to select, for each element, the material point which is the most critical. This way, the amount of gradient data that is written to the database is reduced, which is important when working with laminate materials.

B2000++ discerns between essential boundary conditions (here: displacement constraints) and natural boundary conditions (here: forces). For our model we have to define the essential boundary conditions for (1) creating the required displacement constraints along edge 2 in the x-direction and (2) properly locking the structure along the edges.

Essential boundary conditions ebc along patch edges, i.e displacement constraints, can be added by referencing the proper epatch element (here: edge e2). MDL automatically generates the mesh node and element lists. They are stored on the model database and can be referenced by ebc and nbc boundary condition definitions.

The first ebc block 1 locks the structure. Since the edges are free to rotate we only lock the displacement DOFs, i.e DOFs UX, UY, or UZ::

ebc 1
   # Zero-value constraints
   value  0. dof [  UY UZ] epatch 1 e1  epatch 1 e2  epatch 1 e3   epatch 1 e4
   value  0. dof [UX     ] epatch 1 e4
end

The first line of the ebc block locks the y- and z-displacement DOFs of all edge nodes. The second line locks the x-displacement DOFs of the edge opposite to the moving edge (e2).

The second ebc block 2 applies a value for the x-displacement of all nodes of the edge e2:

ebc 2
   # Forced displacement
   value -1. dof  UX epatch 1 e2
end

The links to the ebc sets will be added to the case block, as will be described later.

For the particular configuration, i.e. a flat plate, an ‘imperfection’ has to be defined to trigger the post-buckling behaviour. This can be achieved by defining a small out-ouf-plane force that ‘pushes’ the plate in the desired post-buckling direction (in the present model only one configuration is possible). The out-ouf-plane concentrated force of 1 N is assigned to the shell patch center node (identified with the epatch node p9) with a nbc block:

# Central out-of-plane force to create imperfection
nbc
  set 1
    value 1. dof FZ epatch 1 p9
  end
end

Note that, if node-local coordinates systems are defined (this is typically the case for curved shells), boundary conditions are defined with respect to the local coordinate system of the nodes associated to the DOF. In this example, however, there are no node-local coordinate systems. Hence, all boundary conditions are expressed in the global coordinate system.

All parameters related to an analysis case description are contained in a case block. In out example we will carry out 2 types of analysis, a buckling analysis and a post-buckling analysis.

Buckling Analysis

The case block for the buckling analysis defines the type of analysis, includes the links to the boundary conditions, and, for the buckling analysis type, defines the number of buckling modes to be computed.

case 1
  analysis  linearised_prebuckling
  nmodes    10
  ebc       1
  ebc       2
  gradients 1
end

adir
  case 1
end

In the analysis case description above, the buckling solver is activated for 10 buckling modes (note that selecting only a very small number buckling modes, such as 1 or 2, might lead to numerical problems in the eigensolver). The clamping and the displacement boundary conditions (ebc 1 and ebc 2) are activated.

The buckling analysis is carried out by entering the shell command

prompt> b2000++ buckling.mdl

After the analysis has completed, the b2print_buckling_loads print the critcial loads:

b2print_buckling_loads buckling
Eigenvalues and buckling loads, case 1
Contributing EBC identifiers: 2
Total reaction force F=(-158593, 2.50111e-12, -3.66168e-15), amplitude Famp=158593

        Mode       Lambda   Crit. load
           1      0.01718         2725  *
           2      0.03948         6262  *
           3      0.05028         7973  *
           4      0.06912    1.096e+04  *
           5      0.08614    1.366e+04  *
           6       0.1099    1.743e+04  *
           7       0.1128    1.789e+04  *
           8       0.1335    2.117e+04  *
           9       0.1628    2.582e+04  *
          10       0.1632    2.588e+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

Results can also be viewed with the B2000++ browser. Enter the shell command below and click on “Solutions”, “case 1”, and “Buckling Loads”:

prompt> b2browser buckling

The buckling modes can be plotted with the baspl++ viewer. Enter the following shell command which displays the plate with the first buckling mode:

prompt> baspl++ plot_buckling_mode.py

The following image shows the shapes of some buckling modes:

Buckling modes (amplitude plotted).

Post-Buckling Analysis

The case block for the post-buckling analysis defines the type of analysis, includes the links to the boundary conditions, and, for the buckling analysis type, defines the number of buckling modes to be computed.

case 1
   analysis        nonlinear
   ebc             1
   ebc             2 sfactor 0.1  # up to 0.1mm edge shortening
   nbc             1 sfunction 1. # constant force, not scaled by load factor
   step_size_init  0.01           # to get smoooth pre-buckling path
   step_size_max   0.01           # to get smoooth post-buckling path
   gradients       1              # compute gradients for all increments
   rcfo_restrict epatch 1 e2      # Restrict computation of reaction forces
                                  # to edge e2 nodes
end

adir
  case 1
end

The following comments pertain to the the post-buckling analysis case description above:

• The ebc 2 link includes the sfactor 0.1, which means that the load is incremented up to 0.1 times the ebc 2 values.

• The nbc 1 link includes sfunction 1.0, which means that the load of nbc 1 is applied constant during all increments.

• Setting step_size_init and step_size_max both to 0.01 , i.e. keeping the increments constant, will produce a smooth post-buckling path for the graphs, at the expense of higher computation time.2

• The rcfo_restrict option tells any post-processing tool, such as Simples, to resctrict the computation of reaction forces to the nodes of the list (here: the epatch edge e2 nodes).

Postbuckling analysis is carried out by entering the shell command

prompt> b2000++ postbuckling.mdl

INFO:solver.static_nonlinear:18:22:45.812: Start static nonlinear solver for case 1
INFO:solver:18:22:45.813: Resolving linear constraints: 102 equations.
INFO:solver:18:22:45.813: 0 dependent constraints were eliminated.
INFO:solver:18:22:45.813: Reduction matrix (903 x 1005) has 903 nonzeros.
INFO:solver.static_nonlinear:18:22:45.813: Start increment 1 of stage 1 for time 0 and time step 0.01.
INFO:solver.static_nonlinear:18:22:45.885: Start increment 2 of stage 1 for time 0.01 and time step 0.01.
...
INFO:solver.static_nonlinear:18:22:48.304: Start increment 99 of stage 1 for time 0.98 and time step 0.01.
INFO:solver.static_nonlinear:18:22:48.331: Start increment 100 of stage 1 for time 0.99 and time step 0.01.
INFO:solver.static_nonlinear:18:22:48.363: End static nonlinear analysis for case 1 with success.

After the analysis has completed, the simples script plot_postbuckling_graphs.py

prompt> python plot_postbuckling_graphs.py

plots the graphs below:

Mid-point out-of-plane displacement vs. edge shortening plot.

Total reaction force vs. edge shortening plot. The reaction force is the sum of all reaction forces applied to epatch edge e3.

Maximum \(\epsilon_{xx}\) vs. edge shortening plot.

Solid Mechanics: Dynamic

Cable-Stayed Bridge

Note

Location of example case:

$B2EXAMPLES/dynamic/bridge

The lowest eigenfrequencies and their corresponding eigenmodes of a cable-stayed steel-concrete composite bridge are calculated. A global FE model of the structure, with shell and beam elements, has been established. Although the problem of cable-stayed structures is geometrically non-linear it is approximated here by replacing the cable by beams with a given prestress.

The FE model of the Zaltbommel bridge in the Netherlands has been elaborated by the NLR, the mesh consisting of R.S rods, B2.S.RS beams, and Q4.S.MITC quad shell elements. The rod and beam elements refer to steel, whereas the shell elements refer to concrete.

FE mesh.

The free vibration analysis with B2000++ is controlled with the MDL commands adir and case, requesting the 20 lowest eigenfrequencies:

case 1
  analysis free_vibration
  nmodes 20
  ebc 1
end

adir
  case 1
end

The cases option of adir specifies which cases B2000++ should process (here: process case 1). In its turn, the case 1 description in the cases 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. free_vibration will launch the B2000++ free vibration analysis solver (see also comment below).

• The nmodes option specifies the number of eigenmodes or eigenvalues to be computed. In the absence of a shift the smallest nmodes eigenvalues around 0.0 are computed.

By default B2000++ calculates the real symmetric eigenvalue problem with the Implicitly Restarted Lanczos Iteration variant of the Implicitly Restarted Arnoldi Iteration algorithm implemented in the ARPACK, the ARnoldi PACKage library. In addition, the model mass and modal stiffness for each mode are computed. Results are summarized below, printed with the print_eigenfrequencies.py script. Alternatively, results can be displayed with the b2browser application.

prompt> b2print_eigenfrequencies zaltbommel

Eigenfrequencies, case 1
     Mode         FREQ        OMEGA         MODK         MODM
        1       0.4275        2.686    1.281e+07    1.775e+06
        2       0.4397        2.763    4.674e+07    6.122e+06
        3       0.6407        4.025    2.993e+07    1.847e+06
        4       0.6961        4.374    2.356e+07    1.232e+06
        5        0.701        4.404    2.356e+07    1.214e+06
        6       0.8348        5.245    4.142e+07    1.506e+06
        7       0.9579        6.019    4.239e+07     1.17e+06
        8         1.16        7.288    1.068e+08     2.01e+06
        9        1.217        7.647    6.866e+07    1.174e+06
       10        1.243         7.81    1.648e+08    2.701e+06

Results can be rendered directly with the baspl++ viewer. Eigenmodes are visualized either with baspl++’s built-in modes viewing tool. The viewer script view.py produces the plot shown below:

First fundamental free vibration mode no. 1.

Any other mode i can be rendered with the baspl++ function mode (see below):

mode(i)

view.py also provides a function for producing plots of several modes:

def mode(imode):
    """Render mode 'imode' (int).
    """
    m2.field.ch.mode = imode
    f = m.db["MODE.1.0.0.1.{}".format(imode)].desc["FREQ"][0]
    s.label = 'Zaltbommel bridge, Eigenmode {}, f={:.2f} Hz'.format(imode, f)

def images(i,j):
    """Render modes from `i` to `j` in ``png`` files.
    """
    for k in range(i, j+1):
        mode(k)
        s.export_to_file("mode{:03d}.png".format(k))

Calling the above defined function

images(1,8)

from the baspl++ shell will render modes 1 to and produce 8 png files mode001.png to mode008.png that can be assembled to one single plot with the montage application:

montage *.png -scale 33% -mode Concatenate -tile 2x8 montage.jpg

producing

Modes 1 to 8 assembled in one single plot with the montage application.

Linear Transient Analysis of Beam

Note

Location of example case:

$B2EXAMPLES/dynamic/simple_beam

A simply supported beam is loaded by a central load P applied at time t=0 (see also input file beam.mdl). The model has been chosen such that the response under the given loads is linear and the beam thickness to beam length ratio is small, thus rendering the influence of the shear strain negligible. A theoretical transient solution for this configuration is available, formulated with the Euler beam theory [clough]:

\[Dy(t) = \frac{P_{0}L^{3}}{n^{4}EI} \sum_{n=1,3,5} \frac{1 - \cos(\omega_{i}t)}{i^4}\]

where i is the odd (symmetric around \(y = L/2\)) mode number, E the modulus of elasticity, I the moment of inertia around the bending axis, m the beam mass per length unit, and L the length of the beam, and \(Ft(t) = P_0\).

The force is applied at \(t=0\) and kept constant. The beam as displayed below is modeled with B2.S.RS beam elements, see MDL input file beam.mdl.

Beam model.

Three different analyses will be carried out: A static analysis, a free-vibration analysis, and a transient analysis. Hence, three different analysis cases will be defined. All three cases refer to the same essential boundary conditions (MDL)

ebc 1
  # Eliminate dof's 3,4,5
  value 0.0 dof [UZ RX RY] allnodes
  # Support
  value 0.0 dof [UX UY] epatch 1 P1 dof 2 epatch 1 P2
end

Natural boundary conditions will be applied to the static and transient dynamic analysis:

nbc 1
  value 1000. dof UY node 21
end

The static analysis case 1 includes both the essential and the natural boundary conditions defined above:

case 1
  analysis linear
  nbc 1
  ebc 1
end

A free vibration analysis gives an idea about the frequencies to be expected, at the same time checking the free vibration analysis solution and preparing the input for the modal analysis. Theoretical free vibration modes and their corresponding circular frequencies \(\omega\) are available with

\[\omega_{i} = i^{2}\pi^{2}\frac{\sqrt{EI}}{mL^{4}}\]

The free vibration analysis case is defined by specifying the essential boundary condition set and the number of (lowest) eigenvalues to be computed:

case 2
  analysis free_vibration
  nmodes 10
  ebc 1
end

The free vibration analysis solution is obtained by explicitly requiring case 2 in the adir command:

adir
  case 2
end

Eigenvalues are obtained by executing b2000++ and by printing them with the Simples print_eigenfrequencies.py script:

$ b2print_eigenfrequencies --case=2 beam

Eigenfrequencies, case 2
Mode    Frequency        Omega      Modal_K      Modal_M
   1         10.5        65.96        625.2     2.72e+06
   2        41.81        262.7        404.6    2.792e+07
   3        93.42          587        184.7    6.363e+07
   4        144.9        910.6        624.8    5.181e+08
   5        164.5         1034        107.7     1.15e+08
   6        254.1         1597        71.99    1.835e+08
   7        361.1         2269        52.58    2.706e+08
   8          435         2733        623.6    4.658e+09
   9        484.2         3043        40.84     3.78e+08
  10        622.6         3912        33.19    5.078e+08

Note that the theoretical circular eigenfrequencies are derived from an Euler beam model. Thus, the shear, which has a certain influence for the present model (although the beam is slender), is not taken into account and explains small differences between theoretical and calculated circular eigenfrequencies.

While the modal solver takes the damping coefficients of each mode \(\zeta\) for the model from the argument list, the transient solver requires the Rayleigh damping coefficients (see [bathe1996] \(\alpha\) and \(\beta\) to be defined:

\[\alpha = \frac{2 \omega_1 \omega_2 (\omega_1 \zeta_2 - \omega_2 \zeta_1)} {\omega_{1}^{2} - \omega_{2}^{2}}\]

\[\beta = \frac{2(\omega_1 \zeta_1 - \omega_2 \zeta_2)}{\omega_{1}^{2} - \omega_{2}^{2}}\]

\(\omega\) and \(\zeta\) are the circular frequencies and the damping percentage, respectively.

The B2000++ transient solver and the modal analysis script msolve are executed without damping or with 10% damping, see the Makefile in the examples directory. The damped transient analysis eventually converges to the static solution, given a damping coefficient >0. The modal analysis test selects the first 3 eigenvalues.

With circular frequencies available from the free vibration analysis, the modal transient analysis can be performed, either with the b2modal program or with the python script msolve (see example directory):

./msolve ${NELEM} ${MIDSTATION} ${T1} ${NSTEPS} 0.0

The time response curve for the beam mid-station y displacement DY (10% damped response) is obtained with the Simples script plot_time_history (see Makefile in the example directory). It produces the svg file dy.svg:

Modal transient analysis: Time-displacement response.

A linear transient analysis without damping or with Rayleigh damping is defined in case3, which is split in the effective case definition (case 3) and the stage definition (case 31). Adding a stage is necessary for specifying the time range (time_step):

case 3
  title 'Transient analysis'
  stage 31 time_step (T)
end

stage 31 integrates from time 0.0 to T, T being a parametrized value (see file beam.mdl).

The stage in its turn defines conditions of the stage:

stage 31
  nbc 1 sfunction '1.'
  ebc 1
  step_size (DT)
end

nbc 1 sfunction '1.'includes the nbc set 1, scaling it to 1.0 for the whole stage. ebc 1 includes the nbc set 1, keeping it constant during the stage (default). multistep_integration_order selects the time integration scheme, here of order 2. step_size defines the time step or increment, (DT) being a parametrized value.

The transient analysis solution is obtained by explicitly requiring case 3 in the adir command:

adir
  case 3
end

The shell command

make all

executes all damped/undamped tests with the direct transient and the modal solvers and produces csv files. These can then be combined to comparison plots. Example:

b2plot_csv --l=center --width=150. --height=50. modal_0.0.csv transient.csv

produces the comparison plot

Compare linear modal and transient curves, mid-station DY displacement, no damping.

and

b2plot_csv --l=center --width=150. --height=50. modal_0.1.csv transient_10.csv

produces the comparison plot

Compare linear modal and transient curves, mid-station DY displacement, 10% damping.

Solid Mechanis: Global-Local

Vibration of Clamped Beam (Craig-Bampton Method)

Location of example case: $B2EXAMPLES/globallocal/beam

TEXT TO BE REVISED STEPS WRONG

The test computes vibration modes of a clamped beam with the Component Mode Synthesis method (also referred to as the Craig-Bampton method) global-local 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 macro-element. Then, the local model uses this reduced system defined by the global model for making a local (detailed) vibration analysis.

The global-local analysis is performed with the run.py script as follows:

1. Condensate the global model global.mdl.

2. Create the input database of the local model local.mdl.

3. 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.

4. Compute vibration modes with the local model.

Eigenfrequencies can now be extracted from the global model:

$ b2print_eigenfrequencies --case=2 global

Eigenfrequencies, case 2
     Mode         FREQ        OMEGA         MODK         MODM
        1        3.726        23.41        480.8       0.8774
        2        11.17        70.21         4326       0.8777
        3         23.4          147         1857       0.0859
        4        65.76        413.2         5381      0.03152
        5        70.07        440.3    1.672e+04      0.08626
        6        129.6        814.2    1.051e+04      0.01586
        7        196.5         1235    4.822e+04      0.03162
        8        216.2         1359    1.728e+04     0.009362
        9        305.9         1922        506.5    0.0001372
       10        326.1         2049    2.491e+04     0.005934

For comparison, eigenfrequencies of the reference model:

$ b2print_eigenfrequencies --case=2 reference

Eigenfrequencies, case 2
     Mode         FREQ        OMEGA         MODK         MODM
        1        3.725        23.41        480.8       0.8775
        2        11.17         70.2         4326       0.8778
        3        23.35        146.7         1858      0.08637
        4        65.38        410.8         5404      0.03203
        5         69.9        439.2    1.673e+04      0.08675
        6        128.1        805.1    1.057e+04       0.0163
        7        195.2         1227    4.873e+04      0.03239
        8        211.9         1331    1.746e+04     0.009849
        9        301.5         1895        524.9    0.0001462
       10        316.7         1990    2.606e+04     0.006582

Modes 2,5 and 7 are displayed in the montage plot below (the modes are extracted and displayed with the baspl++ script view.py and merged with the montage utility.

Some vibration mode shapes. Dots: Local model. Solid line: Global reference model.

I-Profile Beam

An I-profile beam subjected to thermal loading is modelled with solid elements. The 2 flanges and the web are modelled separately and ’tied’ together to form the complete model, referred to as ’tie’ model. Figure 2 displays the shapes of the separate models (B2000++ epatches). It also show the interfaces to be ’tied’, i.e. the lower and upper faces (black dots) of the web (green shape), which are ’tied’ to the corresponding faces of the lower and upper flanges.

Mesh: Element pach contours of lower flange (yellow), web (green), and upper flange(red). Tie interface points shown in black.

To get reference solutions, the same model is assembled to a single continuous mesh with the B2000++ join option. The continuous mesh, however, can be generated only if the nodes of the interface are matching. This is the case for mesh A (see below).

Results are compared to a model with a single matching mesh (join model) and to analytical results from beam theory.

Meshes

Mesh A is selected such that the mesh nodes at the interface are matching, thus allowing for direct comparison of the tie and join models. The figure below displays the continuous mesh A of the structure.

HE (hexahedral) mesh of continuous model (meshA).

Mesh B is selected such that the mesh nodes at the interfaces are not matching and thus can be used with kinematic coupling only. The quality of the solution is then compared to mesh A. The mesh density of the web of mesh B in x-direction is 66% of the one of mesh A.

Coupled Static Heat-Deformation Analysis

The couple static heat and deformation analysis is performed in 2 steps, a heat analysis followed by a deformation analysis retrieving the temperature distribution of the previous heat analysis.

The heat analysis boundary condition sets the upper flange of the structure to 100 degC and the lower on to 0 degC.

Several deformation analysis boundary conditions are studied:

• Condition C1 clamps the beam at x=0. This creates bending with the give temperature distribution. It is interesting to note that (1) the deformations for this case, at least with beam theory, are independent of the cross section shape, the only intervening parameter being the height and that (2) the beam is stress-free, since it can deform freely. This case has a simple beam theory solution for the beam tip displacement.

• Condition C2 clamps the beam at x=0 and X=L, thus generating stresses.

Results for condition C1 are summarized in the table below, demonstrating the efficiency of the kinematic coupling method, the tip displacements obtained with the 3 methods being practically identical.

Deformation solutions condition C1: Tip displacement in z direction, at y=0.

Method	z-displacement
Mesh A: Join (matching)	-0.0128
Mesh A: Tie (kinematic coupling, matching meshes)	-0.0127
Mesh B: Tie (kinematic coupling, non matching meshes)	-0.0124
Beam theory	-0.0115

Global solutions for the heat analysis and the subsequent deformation analysis for condition C1 and mesh B non-matching) are displayed in the figures below. Stresses in the whole model are 0, as expected.

Heat analysis solution (temperatures) of mesh B.

Deformation solution condition C1 rof mesh B.

Deformation solution condition C1 for mesh B details showing detail of non-matching grid between flange and web.

Solid Mechanis: Static

Cable Truss

Location of example case: $B2EXAMPLES/static/cable_truss

A truss composed of prestressed cables deforms under its own weight (simulated here by forces concentrated at the truss free nodes) and under an additional single force applied at one node. The model has been presented by Kar et.al. [kar] . The non-linear analysis can be executed in one or several stages being executed one after the other in the order they are enumerated in the case definition. Each stage references a case specifying sets of boundary conditions and analysis parameters and the solution will be incremented to a load or time parameter of 1.0 unless specified otherwise by the sfactor option of the stage.

Cable truss mesh. Nodes are marked branch:node, where branch is the branch identifier and node the external node identifier. Boundary conditions (red) are marked LLL, i.e displacements locked in all 3 directions at nodes 1, 2, 3, 6, 7, 10, 11, 12).

Three different cases are presented, all three leading to the same result. The problem must be solved in 4 steps (see also note below):

1. Run the input processor b2ip++, defining all 3 cases (see input file cable.mdl.

2. Run for case 1 with the solver option adir cases=1.

3. Run for case 2 with the solver option adir cases=2.

4. Run for case 3 with the solver option adir cases=3.

The database now contains the 3 runs and the load-stress curves can be extracted (see diagrams below).

Note

The non-linear B2000++ solver can treat more than one nonlinear case in one single run. However, the cases are executed one by one and the load factors (times) are added up. This not what is wanted in this examples: All 3 cases must run up to load factor (time) 1.0! Therefore the problem must be solved in 4 steps:

Case 1: The analysis is performed in one stage, with all boundary conditions (dead weight plus load) applied simultaneously. Case 1 tries to solve for all boundary conditions in one stage, with all natural boundary conditions (dead weight plus load) contained in one set: Essential boundary conditions are defined in set 1 (option ebc 1), all natural boundary conditions are defined in natural boundary condition set 3 (option nbc 3), and the solver will try to solve for the final load multiplier value of the stage, i.e. (default) 1.0 with step_size 1.0. In case the step size is too large, steps can become smaller, but not less that 0.01 (option step_size_min 0.01). :

# External forces applied (dead weight and load)

nbc 3 # All forces applied simultaneously
   dof 3 value 100. nodes 9
   value 8. nodes [4 5 8 9]
end

# Case definition

case 1
   title                    'Dead weight and load together in 1 set'
   analysis                 nonlinear
   increment_control_type   load
   ebc                      1
   nbc                      3
   step_size_init 0.1
   step_size_min 0.01
   step_size_max 0.1
end

# Analysis directive
adir
   case 1
end

The output from B2000++ is listed below (it is obtained with the B2000+ command option -l):

b2000++ -l 'info of all in err' -adir 'cases=1' b2test.b2m
INFO:all:10:45:58.554: Start
INFO:solver.static_nonlinear:10:45:58.555: Start static nonlinear solver for case 1
INFO:domain:10:45:58.555: Total number of DOfs: 36.
INFO:solver.static_nonlinear:10:45:58.556: Start prediction of increment 1 of stage id 1 for time 0.
INFO:solver.static_nonlinear:10:45:58.591: Start correction of increment 1 of stage id 1 for time 1.
INFO:solver.static_nonlinear:10:45:58.595: End static nonlinear analysis for case 1 with success.
INFO:all:10:45:58.595: End of execution

Case 2: Again, the analysis is performed in one stage, with natural boundary conditions (dead weight plus load) applied simultaneously, but enumerated separately:

# Boundary conditions

nbc 1 # Dead weight
  dof FZ value 8. nodes [4 5 8 9]
end

nbc 2 # Applied load
  dof FZ value 100. nodes 9
end

# Case definition

case 2
  title                  'Dead weight and load together in 2 sets'
  analysis               nonlinear
  ebc                    1
  nbc                    1 # Own weight
  nbc                    2 # Single force
   step_size_init 0.1
   step_size_min 0.01
   step_size_max 0.1
end

# Analysis directive

adir
  case 2
end

Case 3: The analysis is performed in two stages, one for the dead weight (case 31) and one for the load (case 32):

# Boundary conditions

nbc
   set 1 # Dead weight
       dof 3 value 8. nodes [4 5 8 9]
   end
   set 2 # Single force
       dof 3 value 100. nodes 9
   end
end

# Case definition

case 3
  title                   'Stage 1 (dead weight) and stage 2 (load)'
  analysis                nonlinear
  stage                  31
  stage                  32
end

stage 31
  title                   'Stage 1 (own weight) and stage 2 (load)'
  ebc                     1
  nbc                     1
   step_size_init 0.1
   step_size_min 0.01
   step_size_max 0.1
end

stage 32
  title                   'Stage 2 (load)'
  ebc                     1
  nbc                     2
   step_size_init 0.1
   step_size_min 0.01
   step_size_max 0.1
end

# Analysis directive

adir
  case 3
end

Case 3 solves the stages 31 and 31 in the order they are enumerated in the case block. Again, essential boundary conditions are defined in set 1 (option ebc 1) which has to specified for each stage! The stage with the own weight, i.e. case 31, refers to the natural boundary condition set 31 (option nbc 31), which is solved up to the load factor 1.0. The solver then continues with the next stage, i.e. case 32, which refers to the natural boundary condition set 32 (option nbc 32), solving for the final load multiplier value of that stage, i.e. (default) 1.0.

All 3 cases give the same result for the final load factor, although the path to the solution is not same - it depends on the load stages. Note that case 3 goes to a final load factor of 2 because of the 2 stages, while cases 1 and 2 go to the final load factor of 1.0.

Displacement of node 9 and stresses in element 10 for load case 1.

Displacement of node 9 and stresses in element 10 for load case 2.

Displacement of node 9 and stresses in element 10 for load case 3.

Deformations can be viewed with baspl++, the viewer script view.py producing the deformation plot with the cable stresses shown below (case 1). To run the viewer from the shell, type

baspl++ -t view.py

Deformed cable truss shape and cable stress (colored elements) and applied force.

In addition, node displacements of all cases can be printed with the B2000++ model browser b2browser.

b2browser cable

and click Solutions, Case 1, Displacements, Cycle 10 (the final cycle):

Displacements and Rotations (DISP), dataset "DISP.1.10.0.1", cycle=10, subcycle=0, case=1

      EID        UX           UY           UZ           RX           RY           RZ       Amplitude
        1 G           0            0            0            0
        2 G           0            0            0            0
        3 G           0            0            0            0
        4 G     -0.3509      -0.3509        0.887        1.016
        5 G     -0.5581        3.269       -2.816        4.351
        6 G           0            0            0            0
        7 G           0            0            0            0
        8 G       3.269      -0.5581       -2.816        4.351
        9 G       3.615        3.615        16.31        17.09
       10 G           0            0            0            0
       11 G           0            0            0            0
       12 G           0            0            0            0
Largest amplitude=17.0906
Absolute printout cutoff value threshold=1e-100

Cylinder with RBE-Elements

Note

Location of example case:

$B2EXAMPLES/static/cylinder_rbe

This example demonstrates RBE elements for defining boundary conditions applied to a section. The model is a cylinder clamped at one end \(z=0\) and loaded at the other end \(z=L\). To introduce the corresponding boundary conditions, RBE elements are introduced:

• The clamping condition is formulated with a RBE element connecting a node at (0,0,0) with all circumferential nodes of the cylinder at \(z=0\). The constraints (ebc) will be applied to that node. The RBE element acts as a rigid diaphragm in the xy-plane.

• The load condition is formulated with a RBE element connecting a node at (0,0,L) with all circumferential nodes of the cylinder at \(z=L\). The force (nbc) will be applied to that node. The RBE element acts as a rigid diaphragm in the xy-plane.

Cylinder FE mesh: Q9 shell element mesh (blue) and RBE mesh (red).

Dimensions and material parameters of the cylinder model:

L = 1.0    # Length of cylinder [m]
R = 0.3    # Radius of cylinder [m]
T = 0.0005 # Thickness of cylinder [m]
Fx = 1.e5  # Axial load [N]
E = 72.4e9 # Young's modulus [N/m2]

Cylinder model: Deformation amplitude, traction case.

Cylinder model: Deformation amplitude of RBE element, traction case.

One Element Laminate Material Example

Location of example case: $B2EXAMPLES/static/laminate1

This one-element model consists of a rectangular four-node Q4 shell element in the global x-y plane The model is “statically determinate”, i.e the number of constrained degrees of freedom is equal to the number of rigid body movements, i.e 3 translations and 3 rotations. The plate is pulled in the element x-direction (fiber 1 direction, see figure below).

The purpose of this example is to

• Demonstrate laminate material in shells with a one layer material with 0°, 90° fiber orientations, and with a four layer 0°/90°/90°/0° material.

• Explain the basic concepts of coordinate systems.

Forces (natural boundary conditions) are added as point forces in order to illustrate the concept of node-local coordinate systems. The total force applied is selected such that the stress in fiber direction 1 becomes 500 MPa:

\[F_{tot}=\sigma_{11}*B*T = 500*2.3*0.15 = 194 N\]

The first test 0.mdl defines a model in the branch-global (=global-global) x-y plane:

• The branch-global system is identical to the global-global system (no branch orientation specified).

• The coordinates are formulated in the branch-global system.

• The element local x-axis coincides with the branch-global x-axis.

• The node-local directions of the degrees-of-freedom coincide with the branch-global direction.

• The zero constraints (essential boundary conditions, ebc) do not need node-local transformations.

• The applied forces \(F=0.5*F_{tot}\) in the branch-global x-direction do not need node-local transformations.

One-element laminate test: Element local x-axis coincides with branch-global x-axis. Gray: mesh, blue: deformed mesh, red: constrains, green: forces.

The second test 30.mdl defines a model in the branch-global (=global-global) x-y plane with the element local x-axis rotated by 30° around the branch-global z-axis with respect to the branch-global x-axis:

• The branch-global system is identical to the global-global system (no branch orientation specified).

• The coordinates are formulated in the branch-global system.

• The orientation of the values at the mesh nodes (forces, constraint) are formulated in a node-local orientation system with the MDL transformation command

  transformations
      1 cartesian 0 0 0  0 0 1  (cos30) (sin30) 0
  end
  

  and referred to by the dof_ref parameter of the MDL nodes command:

  nodes
      dof_ref 1 # refers to transformation 1
      1 0.0 0.0 0.0
      2 7.53442101 4.35 0.0
      3 6.38442101 6.34185843 0.0
      4 -1.15 1.99185843 0.0
  end
  

• The zero constraints (essential boundary conditions, ebc) are formulated in the node-local system, because the nodes involved have a local transformation.

• The applied forces \(F=0.5*F_{tot}\) are formulated in the node-local system, because the nodes involved have a local transformation.

One-element laminate test: Model rotated 30° around the branch-global z-axis. Gray: mesh, blue: deformed mesh, red: constrains, green: forces.

The third test 30btrf.mdl defines a model in a branch-global

transformed with respect to the global-global system by a branch_orientation command:

• The branch-global system is rotated with the MDL branch_orientation command by 30° around the global-global z-axis:

  branch_orientation
      rotate axis z angle 30
  end
  

• The coordinates are formulated in the branch-global system and are identical to the coordinates of model 0.mdl.

• The element local x-axis coincides with the branch-global x-axis.

• The node-local directions of the degrees-of-freedom coincide with the branch-global direction (no node-local transformation).

• The zero constraints (essential boundary conditions, ebc) do not need node-local transformations.

• The applied forces \(F=0.5*F_{tot}\) in the branch-global x-direction do not need node-local transformations.

Comparison of the results of the three models (the solutions of the analysis) is discussed. To resume, please note that:

• DOF fields, i.e. displacements and rotations, are stored on the B2000++ model database with respect to the branch-global coordinate system. Values at nodes with node-local coordinate systems (orientations) are transformed to the branch-global orientation.

• Element stress and strain tensors as well as any other sampling point field values, i.e element-related derived quantities, are stored on the B2000++ model database with respect to the branch-global coordinate system.

Displacements of nodes 2 and 3 as stored on the B2000++ model database are compared in the table below:

Comparison of displacements

   Model  EID           DX           DY    Amplitude
     "0"    2      0.02962            0      0.02962
            3      0.02962    -0.002349      0.02971

    "30"    2      0.02565      0.01481      0.02962
            3      0.02683      0.01278      0.02971

"30btrf"    2      0.02962            0      0.02962
            3      0.02962    -0.002349      0.02971

Displacements of model 0 and 30btrf are identical, because both models are formulated in the same branch-local configuration.

Stresses at specific points can be plotted or printed with the b2plot_stress_stack script.

script. Example: Plot stresses through the shell thickness z' at integration point 1 of element 1.

b2plot_stress_stack -e 1 -p 2 --case 1 --print -f STRESS 30;

Model 30: Stresses through shell thickness z' at integration point 1 of element 1.

Print stresses through the shell thickness z' at integration point 1 of element 1.

b2plot_stress_stack -e 1 -p 2 --case 1 -f STRESS --print 30

Layer          z   Sigma_xx   Sigma_yy   Sigma_zz   Sigma_xy   Sigma_yz   Sigma_xz
    0      -0.08        701        246   8.98e-35        394   3.56e-17  -2.06e-17
    0      -0.04        701        246   2.25e-35        394   1.78e-17  -1.03e-17
    1      -0.02       48.7        3.9  -3.01e-36       38.8  -6.57e-18   3.79e-18
    2       0.02       48.7        3.9  -3.01e-36       38.8   6.57e-18  -3.79e-18
    3       0.04        701        246   2.25e-35        394  -1.78e-17   1.03e-17
    3       0.08        701        246   8.98e-35        394  -3.56e-17   2.06e-17

Scordelis-Lo Roof Linear Analysis

Location of example case: $B2EXAMPLES/static/scordelis_lo_roof

The Scordelis-Lo roof is a classical test case for linear static analysis named after the authors who first reported it [scordelis]. The test has been proposed as a standard FE test by MacNeal and Harder [macneal1985] and it can be found in the NAFEMS test case description.

The structure consist of a concrete cylindrical shell roof supported by rigid diaphragms on both curved sides (see figure below), and it is loaded by a surface traction of 90psf in the negative y-direction. The total length of the cylinder is 50ft, the radius of the cylinder 25ft, the half opening angle 40°, and the thickness 0.25 ft. The material properties are isotropic, with the elastic modulus of 4.32e8psf and the Poisson number set to 0 or, in some references, to 0.17. Von Mises failure stress is 300000psf.

Scordelis-Lo Roof: Mesh.

The diaphragms are defined such that a simply supported boundary condition results: Displacements in x- and y-directions are suppressed, the suppression of displacements in the z-direction being assumed by the symmetry condition along \(P_4\)-\(P_1\).

The test compares the linear static solution, specifically the displacement in the y-direction of point \(P_4\), against known solutions (the reference solution for point \(P_4\) is -0.3024).

The solution of this problem with B2000++ consist of the following steps:

• Create the finite element model with the MDL language describing the problem in a text file.

• Solve the problem by submitting the text file to the B2000++ solver. The solver creates a B2000++ database.

• Examine the results by reading the B2000++ database.

A finite element model for B2000++ is created by describing the problem geometry, creating a finite element mesh, adding boundary conditions, and adding solver instructions. The complete model can be formulated with the B2000++ Model Description Language (MDL). Alternatively, parts of the model, such as the mesh and the boundary conditions, can be produced with an external program and brought in the B2000++ input file.

The Scordelis-Lo roof problem can easily be produced with the B2000++ Model Description Language. The necessary steps will now be explained. All steps result in instructions added with an editor to a single B2000++ input file.

Geometry and mesh: The cylinder geometry and the mesh produced by the geometry can be described with an epatch (element patch). Element patches allow for generating simple building blocks, such as plates, cylinders, etc. The Scordelis-Lo cylindrical roof is described by an epatch of type cylinder, which implies the creation of a surface and a surface mesh:

epatch 1
    geometry cylinder
    type (ELTYPE)     # Element typ
    ne1 (NE) ne2 (NE) # N. of elements along sides
    thickness 0.25    # Thickness of shell
    mid 1             # Material identifier
    phi1 50.0         # Opening angle of half cylinder
    phi2 90.0         # Closing angle of half cylinder
    radius 25.0       # Radius of cylinder
    length 25         # Half length of cylinder
    # local edges means that all DOF's of the nodes along the edges are
    # expressed in local, i.e. cylindrical, coordinates.
    if (local == 1) {local edges}
end

Note that the epatch must contain the element type (eltype keyword), the number of elements to be generated along each side of the patch (ne1 and ne2 keywords), and material reference (mid keyword). The definition of the epatch is further illustrated in the figure below, with the position of the vertices P1 to P4 and the edges E1 to E4.

In this particular configuration the local and global orientations are the same, i.e. no need to define cylindrical coordinates, unless the resulting displacements, rotations, forces, and moments should be expressed in cylindrical coordinates (see instruction if (local == 1) {local edges} above).

Scordelis-Lo Roof: epatch vertex points and edge definitions. See also the generated mesh.

Material: The material is the same as the one in the original problem description. The definition is straightforward with the material'' command, which identifies the material with id ``1

material 1
    type isotropic
    e 4.32e8       # Young's modulus
    nu 0.          # Poisson number
    density 90
    failure von_mises
      r 300000
    end
end

Boundary conditions: Since all boundary conditions apply either to edges of the epatch or to the epatch surface they can be generated for all nodes of an edge by referring to the epatch edges which in their turn are automatically generated. Essential boundary conditions (here: displacement constraints) are defined with the ebc command:

ebc 1
    value 0 # Values to be assigned to DOF's
    # Diaphragm = Supported edge (global=local orientation) along
    # edge e1 from P1 to P2
    dof [UX UY] epatch 1 e1
    # Symmetry condition along edge e2 from P2 to P3
    dof [UX RY RZ] epatch 1 e2
    # Symmetry condition along edge e3 from P3 to P4
    dof [UZ RX RY] epatch 1 e3
end

The figure below shows the essential boundary conditions as defined above: At each node of the edges a pattern XXXXXXX is plotted, where the position of X in the pattern designates the degree of freedom (DOF) to be affected and X taking the values F (free, i.e. unconstrained), L (locked, i.e. constrained to 0), or V (constrained to a value).

Scordelis-Lo Roof: Essential boundary conditions generated by ebc 1.

The surface pressure loads are applied with the natural boundary condition command nbc. Note that the reference frame for the surface tractions is the branch reference frame and thus has to be specified explicitly with the system option, because the default reference frame is the local reference frame relative to the element surface.

nbc 1 type surface_tractions
    system branch # Branch referential (default is local!)
    surface_tractions 0 -90 0
    epatch 1 f7   # Apply to mid-surface of shell!
end

What remains now to be specified are (1) definition of the case, comprising the type of analysis (2) the boundary condition sets to be applied to the case (3) the analysis directives, i.e. the cases to be solved for:

case 1
    analysis linear
    nbc 1
    ebc 1
end

adir
    case 1
end

The analysis is run with the shell command:

b2000++ scordelis_lo_roof.mdl

which generates the B2000++ database scordelis_lo_roof.b2m. Specifying the log option to the b2000++ command prints a short summary of the steps which the solver performs:

$ b2000++ -l info scordelis_lo_roof.mdl
INFO:all:13:57:24.656: Start
INFO:solver.static_linear:13:57:24.658: Start static linear solver
INFO:domain:13:57:24.658: Total number of DOfs: 422.
INFO:solver.static_linear:13:57:24.659: Element matrix assembly
INFO:solver.static_linear:13:57:24.660: Resolve the linear problem
INFO:solver.static_linear:13:57:24.664: Compute gradients and reaction forces
INFO:solver.static_linear:13:57:24.664: End of static linear solver
INFO:all:13:57:24.664: End of execution

Results

The summary of the analysis results for the 8 by 8 element mesh is displayed in the table below for element type Q4.S.MITC.E4:

Q4.S.MITC.E4 results (8 by 8 element mesh)

What	Reported	B2000++	Error (%)
Displacement DY at point P3 [ft]	-0.3024	-0.304196	0.6

Results can be examined with several B2000++ tools, such as the baspl++ post-processor, the Simples scripting tools, the *b2print_ text output tools, or the B2000++ data browser. We start with the baspl++ viewer to get an image of deformed shape of the structure. The figure below is obtained by launching baspl++ from the shell with the viewer script view.py:

baspl++ -t plot_displacements.py

and shows the deformed structure (colored) and the initial undeformed one (black dashed grid). Please read the file view.py.

Scordelis-Lo Roof: Displacement field.

baspl++ -t plot_stresses.py

plots the per element maximum von Mises comparison stress:

Scordelis-Lo Roof: von Mises stresses (max of element).

baspl++ -t plot_stresses_sampling.py

plots the von Mises comparison stress at each integration (sampling) point in the form of spheres, where color and sphere size are proportional to the intensity:

Scordelis-Lo Roof: von Mises stresses (all sampling points).

Values can be extracted with Simples scripts. The script print_displacement.py extracts the displacement at point P3 and prints it:

python print_displacement.py

producing the output

Displacement at point P3 (node 273): DY=-0.301733
Reference solution DYREF=-0.3024
Error=-0.22%

and the script print_stresses.py extracts and prints stress values at selected elements prints it:

python3 print_vmises.py

Max failure index fmax=1.124541, element=241

Data can be viewed with the B2000++ data browser. Example: View the failure index:

Scordelis-Lo roof: View failure index values with the B2000++ browser.

A convergence study with Simples scripting is very easily performed. Check the scripts convergence.py and convergence_all.py. Results of the convergence study with different element types is displayed below:

Scordelis-Lo roof: Convergence graph for selected element types.

Linear Analysis of Plate with Subcases

Location of example case: $B2EXAMPLES/static/subcases

Th example illustrates the concept of subcases (see case definition). During linear analysis, subcases can be applied to a common set of essential boundary conditions, thus reducing solution time, because the constraints are applied once and the global matrix is factored once. Each subcase is a new right-hand side and needs back-substitution only.

The model is the same as the Buckling of Composite Plate example: A simply-supported thin quadratic laminated plate of 500 mm by 500 mm is loaded with several load cases and several essential boundary conditions.

The first load case (case 1) constrains the plate to simply supported condition (essential boundary condition 1). Under this condition, 2 natural boundary condition (nbc) orce sets are applied as subcases:

case 1
    analysis linear
    gradients 1
    ebc 1
    subcase 1 # Subcase 1: Pressure load (nbc 1)
        nbc 1
    end
    subcase 2 # Point load applied at the plate's mid-point (nbc 2)
        nbc 2
    end
end

The second load case constrains the plate to clamped along edge 4. Under this condition, 2 natural boundary condition sets (forces) are applied as subcases:

• Subcase 1: Pressure load (nbc 1).

case 1
    analysis linear
    gradients 1
    ebc 1
    subcase 1 #Subcase 1: Pressure load (nbc 1)
        nbc 1
    end
    subcase 2 # Two point loads at free vertices creating a torsional
                moment (nbc 3).
        nbc 3
    end
end

Geometry and essential boundary conditions of plate.

Note that B2000++ always solves for all subcases of a case, i.e specifying

adir
    case 1
end

will solve for subcase 1 and subcase 2 of case 1.

The input file plate.mdl contains all definitions for solving the 2 cases with 2 subcases each. The essential and natural boundary condition definitions are not listed here (see MDL input file plate.mdl).

Case 1 defines the 2 subcases solved with the constraint set ebc 1:

case 1
    analysis linear
    ebc 1
    subcase 1
        nbc 1
    end
    subcase 2
        nbc 2
    end
end

Case 2 defines the 2 subcases solved with the constraint set ebc 21:

case 2
    analysis linear
    ebc 2
    subcase 1
        nbc 1
    end
    subcase 2
        nbc 3
    end
end

Finally, adir requests to solve case 1 and case 2:

adir
   cases [1 2]
end

The deformed shapes of the subcases are displayed in the figures below. Note that the selection of the subcase in baspl++ and Simples is obtained with the subcase= or mode= parameter.

Linear Analysis of Plate with Subcases: Amplified deformed shape, case 1, subcase 1.

Linear Analysis of Plate with Subcases: Amplified deformed shape, case 1, subcase 2.

Linear Analysis of Plate with Subcases: Amplified deformed shape, case 2, subcase 1.

Linear Analysis of Plate with Subcases: Amplified deformed shape, case 2, subcase 2.

Heat Transfer: Non-Stationary

Transient Heat Conduction in Soil

Location of example case: $B2EXAMPLES/heat/soil_freezing

The cooling of soil under a constant temperature of the air in the atmosphere and an initial temperature of the soil is calculated for a period of 60 days, the example being taken from the text book by Kreith and Bohn [kreith]. This 3D example demonstrates the step size control of the non-stationary nonlinear heat analysis solver.

The soil is modeled by a row of solid heat conduction elements (HE8 elements) extending 3 m down in the soil (z-direction). The soil is assumed to be semi-infinite, thus, a row of elements is sufficient, with free faces in the x- and y-directions. The face of the model which is in contact with the atmosphere is permanently set to -15 °C, while the rest of the model is set +20 °C for t=0. It is assumed that the conduction properties of the soil are uniform and linear and no convection or radiation between the soil and the atmosphere are taken into account.

The FE mesh is generated with a single epatch, with the face f1 of the epatch being in contact with the air at \(x=0\) (see figure below).

Transient heat conduction in Soil: Mesh and boundary conditions.

The essential boundary conditions are a temperature of -15 on the face which is exposed to the air:

ebc 1
  value -15. dof 1 epatch 1 f1
end

The initial conditions at the mesh nodes specified by means of the dof_init command:

dof_init 1
  value 20 allnodes      # Initially all nodes
  allow_override yes     # Default is 'no'
  value -15. epatch 1 f1 # Overwrite nodes of face to air to match ebc
end

The transient response is calculated with the nonlinear transient solver, and the solution is controlled by the case and stage commands.

By default the solver increments the time parameter from 0.0 to 1.0. In this example the time parameter range is from 0 to 60 days, i.e. 5184000 seconds. The best way of specifying this range is by introducing a case with a stage and by scaling the stage with the time_step option:

case 1
  analysis              dynamic_nonlinear # Although problem is linear
  stage 2 time_step     5184000   # Time is (60 days in seconds)
end

The stage 2 (see listing below) then defines the parameters of the actual case to be computed: The initial conditions and the essential boundary conditions, no natural boundary conditions being defined in this model.

The essential boundary conditions are activated with the ebc option of case 2:

ebc 1 sfunction '1'

By default, non-zero essential and natural boundary conditions are scaled by the analysis time. Specifying sfunction 1 tells the solver to keep the ebc constant in time.

By default, the nonlinear transient uses its error estimator to determine the step size, which results in a good computational efficiency. The option tol_dynamic sets the normalized maximum time integration error. Setting this to a lower value will increase accuracy but will result in a larger number of time steps. Conversely, setting this to a higher value will result in a faster analysis at the expense of accuracy.

stage 2
  dof_init       1               # Initial conditions
  ebc            1 sfunction '1' # Keeps ebcs constant
  tol_dynamic    0.01            # Time integration error control
end

Finally, the adir command specifies the analysis case to be solved:

adir
  case 1
end

The temperature at a given mesh node (at 0.7 m below the surface) as a function of time for several runs with a different dynamic error tolerance is executed with a Simples script (see file run.py):

# Run B2000++ with several time integration error tolerances TOLD,
# extract solutions, and add to list 'solutions' for plotting.

solutions = []
for told in (1000., 1., 0.1):
    os.system('b2000++ -define TOLD=%g soil.mdl' % told)
    model = si.Model('soil')
    x = model.get_steps(key = 'TIME')
    y = model.get_dof_value_of_steps('temperatures', 8, 1)
    model.close()
    solutions.append((x, y, told))

The solution list is then taken up by the graph generation section of run.py, which generates a PNG file:

plotted in the figure below. The dynamic error tolerance controls the time step: When the dynamic error tolerance is high (red markers in the graph below), less steps are performed, at the expense of accuracy. Note that the graph below is created with the Simples tools (see run.py script).

Transient Heat Conduction in Soil: Temperature T(t) at z=0.7[m]. The dynamic error tolerance controls the time step: When the dynamic error tolerance is high (red markers), less steps are performed, at the expense of accuracy.

Note: A constant step size - at the expense of accuracy - can be achieved by setting the step_size_init, step_size_max, and step_size_min options to the same fixed value, and by specifying a high tolerance for the time integration error with the tol_dynamic option:

stage 2
  step_size_init 86400           # 1 day
  step_size_max  86400           # Keeps step size constant
  step_size_min  86400           # Keeps step size constant
  tol_dynamic    10000           # Time integration error control
  dof_init       1               # Initial conditions
  ebc            1 sfunction '1' # Keeps ebcs constant
end

The analytical function can be found in the textbook by Kreith and Bohn [kreith] and it is programmed in the run.py script.