# 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): 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. 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')  ### 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.

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.

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

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

### 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. 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: 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: ## 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.

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

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:

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

### 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. 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:

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

and

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


produces the comparison plot

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

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
analysis                 nonlinear
ebc                      1
nbc                      3
step_size_init 0.1
step_size_min 0.01
step_size_max 0.1
end

# Analysis directive
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

dof FZ value 8. nodes [4 5 8 9]
end

dof FZ value 100. nodes 9
end

# Case definition

case 2
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

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
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
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
ebc                     1
nbc                     2
step_size_init 0.1
step_size_min 0.01
step_size_max 0.1
end

# Analysis directive

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.

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


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. 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]  ### 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.

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.

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;


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

baspl++ -t plot_stresses.py


plots the per element maximum von Mises comparison stress:

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:

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:

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:

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

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

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.