# Selected Topics

## Analysis Cases

The Analysis Cases case refer to the ingredients to be used for an analysis. Specifically, they enumerate the boundary condition sets and other constraint sets as well as any additional sets required for a specific analysis. In addition, the analysis type and the corresponding analysis parameters, as well as numerical parameters, may be specified (these will override the analysis type and parameters defined in the adir analysis directives).

The case mechanism work as follows: The keyword `cases`

of the
ADIR dataset contains the list of case
identifiers to be processed or solved for. This list is subset of the
defined cases.

Each case is defined by the MDL case command and
described on the database under the identifier CASE.id, where `id`

is a positive integer.

**Example** of a case definition

A case with identifier `7`

is defined by

A reference to an essential boundary condition set

`3`

.A reference to 2 natural boundary condition sets

`15`

and`7`

A reference to a linear constraint set

A reference to a branch connectivity set.

The solver will then produce (1) the constraint set for case 7 with the EBC dataset set, the linear constraint set LINC dataset, and the JCTS dataset branch connectivity set and (2) the right hand side (natural boundary condition) with the sum of the NBC dataset sets 15 and 7.

## Solver Algorithms

### Linear Static Solver Numerical Algorithm

The static linear problem exposed by the model abstraction layer of
*B2000++* is:

\[\begin{split}\left\{ \begin{array}{l} {F_{int}u + f_{int} = F_{ext}u + f_{ext}} \\ {u = {Rv} + r} \\ {Lu + l = 0} \\ \end{array} \right.\end{split}\]

where \(F_{int}u + f_{int}\) is the linearized internal value of the domain (the internal force in a displacement based formulation), \(F_{ext}u + f_{ext}\) is the linearized natural boundary condition (the external force in a displacement based formulation), \(R\) and \(r\) is the linearized model reduction, \(L\) and \(l\) is the linearized essential boundary condition.

Using the model reduction equation, the problem becomes:

Using the Lagrange Multiplier Adjunction method for imposing the essential boundary conditions, the problem becomes:

This linear system is solved and the solution \(u\) is then obtained form the linearized model reduction equation.

The difference of the natural boundary condition with the internal value of the domain (this is the reaction force in a displacement based formulation) is then:

### Nonlinear Solver Numerical Algorithm

A incremental-iterative algorithm, also called predictor-corrector method, is used to solve this problem. A stage is broken down into increments and for each increment a corrector phase is applied.

The static nonlinear problem exposed by the model abstraction layer of
*B2000++* is:

After imposing the reduction, the problem becomes:

The potential energy of the unconstrained problem noted \(\Pi\left( {u,\lambda} \right)\) is defined by this first-order variation.

Thus we have:

To impose the constraint, define Lagrange multipliers \(\Lambda\) to form the Lagrangian:

The first-order derivatives are:

where \(c\left( v,\Lambda,\lambda \right)\) is the increment constraint condition that defines the increment control strategy.

The second-order derivatives are:

### Nonlinear Dynamic Solver Numerical Algorithm

The dynamic nonlinear problem is:

The general form of a LMS (linear multi-step) scheme for a normalized first-order ODE (ordinary differential equation) of the form:

is:

where

and

Adapted to a second-order ODE (ordinary-differential equation) on u where:

we obtain:

where

At each time integration step, the problem to solve is:

The first-derivative of the potential energy of the unconstrained problem is:

To impose the constraint, add Lagrange multipliers to form the Lagrangian:

The first-order derivatives are:

The second-order derivatives are:

### Linearised Pre-Buckling Solver Numerical Algorithm

The static nonlinear problem exposed by the model abstraction layer of
*B2000++* is:

A linearized static solution is first computed.

### Free Vibration Solver Numerical Algorithm

The undamped free vibration problem exposed by the model abstraction
layer of *B2000++* is:

Where \(Ku\) is the linearized elastic force vector, \(M\overset{\operatorname{˙˙}}{u}\) is the linearized inertial force vector, \(R\) is the linearized reduction matrix and \(L\) is the linearized essential boundary condition matrix.

Using the model reduction equation, the problem becomes:

Using the Lagrange Multiplier Adjunction method for imposing the essential boundary conditions, the problem becomes:

The free vibration frequencies and modes are then found by solving the following eigenvalue problem:

This eigenvalue problem is then solved using the Implicitly Restarted Arnoldi Iteration algorithm [Sorensen92] (or the Implicitly Restarted Lanczos Iteration variant if the problem is real and symmetric) implemented in the ARPACK library.

The Implicitly Restarted Arnoldi Iteration algorithm computes the higher eigenvalues of the problem. To compute the lower eigenvalues or the eigenvalues near a given position \(\sigma\), a spectral shift value is applied to the eigenvalue problem that becomes:

The eigenvalues can then by obtained using the relation:

Finally, the eigenfrequencies are:

and the eigenmodes are:

### Frequency-Dependent Vibration Solver Numerical Algorithm

The frequency-dependent undamped free vibration problem exposed by the
model abstraction layer of *B2000++* is:

Where \(K\left( \lambda \right)u\) is the frequency dependent linearized elastic force vector, \(M\left( \lambda \right)\overset{\operatorname{˙˙}}{u}\) is the frequency dependent linearized inertial force vector, \(R\) is the linearized reduction matrix and \(L\) is the linearized essential boundary condition matrix.

**First stage:.**

All eigenfrequencies are computed for a given frequency using the same algorithm than the one used for the Undamped Free Vibration Solver. These approximated eigenfrequencies will then be used in the second stage to compute the exact eigenfrequencies.

**Second stage.**

For each of the approximated eigenfrequencies, the following iterative algorithm is applied.

Computation of \(K\left( \lambda \right)\) and \(M\left( \lambda \right)\) with the last computed approximated eigenfrequency value.

Compute the eigenfrequency value near the approximated eigenfrequency by solving the standard eigenvalue problem using the same algorithm than the one used for the Frequency-Dependent Undamped Free-Vibration Solver (B2000++ Pro).

If the iterative eigenvalue solver does not converge, the number of Lanczos vectors used is increased until this algorithm converges.

If the new computed eigenfrequency is not nearly equal to the approximated eigenfrequency (\(\frac{\lambda_{i} - \lambda_{i - 1}}{\lambda_{i}} > \text{tol}\)), restart the iteration at step 1 with the approximated eigenfrequency just calculated.

If the exact eigenfrequency value is too far from the approximated eigenfrequency (\(\frac{\lambda_{i} - \lambda_{0}}{\lambda_{i}} < \text{tol\_predictor}\)), the list of all the eigenfrequencies is recomputed with the current eigenfrequency.

If the new approximated eigenfrequency does not correspond to the eigenfrequency computed in step 3, the eigenfrequency of the current mode is recomputed. Go to step 1.

Store the eigenfrequency and the associated mode in the database and go to the next eigenfrequency computation at step 1.

This numerical method only works if the eigenfrequencies of the problem are well separated.

## Artificial Damping

The MDL command `residue_function_type artificial_damping`

specifies
that artificial damping should be performed. The artificial damping
method performs viscous damping. The “velocities” are defined by the
ratio of the last converged step’s incremental displacements and the
size of the last converged step’s load increment.

The global viscous damping matrix can be constructed in two ways. By default, the global viscosity matrix C is computed from the global mass matrix M:

The coefficient \(\alpha\) is determined in the first load step,
and such that the ratio of dissipated energy (due to the artificial
damping) and of the total strain energy will be equal to the value of
the parameter `dissipated_energy_fraction`

(default 1.0e-4). During
all subsequent load steps, the coefficient \(\alpha\) is kept
constant.

If, on the other hand, one or both of the coefficients
`rayleigh_alpha`

and `rayleigh_beta`

are set, Rayleigh damping
will be performed where the global viscosity matrix C is a linear
combination of the global mass and linear global stiffness matrices M
and K:

Note

The densities of the materials used must not be zero, otherwise the global mass matrix will evaluate to zero. The global mass and global linear stiffness matrices are constant and are evaluated only in the first load step.

The coefficients `rayleigh_alpha`

and `rayleigh_beta`

can be
estimated from the first two (free-)vibration eigenmodes, based on
percentages of critical damping of the associated linear vibration
problem (see also [Bathe96]). Note that, if the circular frequencies
of the two modes are close, the critical damping ratios should be
close as well.

The ``b2print_rayleigh_damping_coefficients`

extracts the
coefficients from a database containing a free-vibration
analysis. Example:

```
b2print_rayleigh_damping_coefficients demo.b2m
```

The same is achieved with the following Python function, which takes the circular frequencies of the first two eigenmodes and the ratios of critical damping, and calculates and prints the coefficient values.

```
def get_alpha_and_beta(omega_1, omega_2, zeta_1=0.001, zeta_2=0.001):
n = (omega_1**2 - omega_2**2)
alpha = 2.0 * omega_1 * omega_2 * (omega_1 * zeta_2 - omega_2 * zeta_1) / n
beta = 2.0 * (omega_1 * zeta_1 - omega_2 * zeta_2) / n
print 'rayleigh_alpha ', alpha
print 'rayleigh_beta ', beta
```

## Shell Elements and the Drill Stiffness

A well-known problem with shell elements is how intersections or junctions (such as the connecting edge between panels and stiffeners) are treated. Pure 5 degrees-of-freedom (DOF) elements are not able to take into account the drill moment and the rotation in the element out-of-plane direction (local z-direction). This violation of the compatibility requirement results in numerical problems that often prevent a successful simulation. Elements with an additional degree-of-freedom are able to maintain compatibility at junctions. However, at nodes where there is no such junction, 6 DOF elements must introduce an artificial drill stiffness term to prevent the second variation from becoming singular.

The MITC elements largely eliminate the need for artificial drill
stiffness by allowing for each node to have either 5 or 6 DOF. The
*B2000++* input processor `b2ip++`

automatically determines which
shell element nodes have 5 DOF, which ones have 6 DOF (see
adir), and which ones of the latter will be
drill-stabilized. This is done such that junction nodes and nodes
where boundary conditions or linear multi-point constraints on one or
several rotational DOF are defined have 6 DOF, while the other nodes
have 5 DOF. Nodes that connect to beam elements also have 6 DOF. The
shell element nodes that are marked to have 6 DOF can be visualized
with the **baspl++** post-processor by selecting the node type
19:

```
if len(sys.argv) != 2:
print 'usage: baspl++ %s database' % sys.argv[0]
sys.exit(1)
m = Model(sys.argv[1])
p1 = NPart(m)
p1.nodes.ntypes = 19
p1.nodes.extract = True
p2 = NPart(m)
p2.edge.show = True
p2.opacity = 0.1
p2.elements.extract = True
```

Nodes whose drill rotation will be stabilized are restricted to the
smallest possible set. For instance, nodes at intersections do not
need stabilization, because the adjacent shell element(s) will provide
the required stiffness.Nodes that are connected to a beam element do
not need stabilization, since the beam element will provide the
required stiffness. For the nodes that will be stabilized, the input
processor b2ip++` will create the node set `b2ip_drills`

which can
be visualized with the following **baspl++** script:

```
if len(sys.argv) != 2:
print 'usage: baspl++ %s database' % sys.argv[0]
sys.exit(1)
m = Model(sys.argv[1])
p1 = NPart(m)
p1.nodes.nodesets = 'b2ip_drills'
p1.nodes.extract = True
p2 = NPart(m)
p2.edge.show = True
p2.opacity = 0.1
p2.elements.extract = True
```

This approach ensures that most of the numerical problems traditionally associated with the 6th degree-of-freedom are resolved. In addition, for the nodes with 5 DOF, node normal vectors are computed by the input processor. Having continuity of the shell normal between elements for proper handling of shell intersections and correctly transferring twisting moments in curved shell models [Hoff95], ensures that adjacent elements are truly compatible along the edges and is a prerequisite for solving problems like the Raasch Challenge (see the Raasch verification test).

This works even when a part of the shell elements is inverted such
that the shell element normal points in the opposite
direction. However, opposite shell element normals may be due to
modeling errors and may cause problems in conjunction with
non-symmetric laminates or traction loads. For this reason, the input
processor `b2ip++`

issues a warning when opposite shell element
normals are detected, and will create the node set `b2ip_opposite`

which can be visualized with the following script:

```
if len(sys.argv) != 2:
print 'usage: baspl++ %s database' % sys.argv[0]
sys.exit(1)
m = Model(sys.argv[1])
p1 = NPart(m)
p1.nodes.nodesets = 'b2ip_opposite'
p1.nodes.extract = True
p2 = NPart(m)
p2.edge.show = True
p2.opacity = 0.1
p2.elements.extract = True
```

## Solutions Fields

When a case is processed, the solver will write the solution data and
meta-information to the database. The *B2000++* solvers produce 3
types of solution fields:

DOF fields: DOF fields are 2-dimensional arrays containing solution data of nodes or elements of the mesh. One single dataset containing the array is defined per branch and case, subcase, cycle, and subcycle. See also B2000++ Database document.

Element fields

Sampling point fields: Element-wise ‘derived’ quantities of the DOF fields. The values are stored element-wise in “Array Tables”. They can contain a variable number of columns for each element, where a column represents a ‘sampling’ or element integration point. In the current version of

*B2000++*the sampling points are collected in groups. One single dataset containing the array is defined per case and cycle.

### Deformation Analysis Solution Fields

For static and transient analysis the solvers of *B2000++* produce the
following solution datasets on the model database:

Displacement and rotation field

`DISP.branch.cycle.0.case`

and, optionally, during transient analysis,`DISPD.branch.cycle.0.case`

, and`DISPDD.branch.cycle.0.case`

.Gradients, such as stresses, strains, beam section forces and section moments.

For linear vibration analysis the solvers of *B2000++* produce the
following solution datasets on the model database

Vibration modes (normalized displacement and rotation fields)

`MODE.branch.cycle.0.case.mode`

. The descriptors of these sets contain the solution attributes, i.e circular frequency, frequency, modal stiffness, modal mass.

For buckling analysis the solvers of *B2000++* produce the following
solution datasets on the model database:

Buckling modes (normalized displacement and rotation fields)

`BMODE.branch.cycle.0.case.mode`

. The descriptors of these sets contain the solution attributes, i.e eigenvalue, critical load.

### Heat-transfer Analysis Solution Fields

The heat-transfer analysis solver of *B2000++* will produce the
following solution datasets on the model database:

Temperature distribution at mesh nodes

`TEMP.branch.cycle.0.case`

Temperature diffusion velocity distribution at mesh nodes

`TEMPD.branch.cycle.0.case`

Heat flux

`HEAT-FLUX.branch.cycle.0.case`

at element sampling pointsHeat capacity

`HEAT-CAPACITY.branch.cycle.0.case`

at element sampling points.

## Coordinate Systems and Transformations

Coordinate systems and related transformations are important for the
understanding of the basic mechanisms of *B2000++*, since the
description of meshes in *B2000++* is substructure (branch)
oriented. A *B2000++* mesh can consist of an arbitrary assemblage of
meshes collected in branches. However, these meshes or branches
themselves may not contain other branches. *B2000++* employs several
Cartesian coordinate systems:

**Global Cartesian coordinate system**

All meshes of all branches are positioned in the global Cartesian coordinate system \(x_g, y_g, z_g\).

**Branch Global Cartesian coordinate system**

The coordinate system \(x_{bg}, y_{bg}, z_{bg}\) in which the mesh
and the solutions of a specific branch are defined and stored on the
*b2000++* model database. If not explicitly specified with the
MDL branch_orientation command, the
branch global coordinate system is identical with the global
coordinate system.

**Node Local Cartesian Base**

The base \(x_{nl}, y_{nl}, z_{nl}\) in which node related
properties are defined stored on the *b2000++* model database, i.e
essential boundary conditions (constraints, see ebc)
and natural boundary conditions (forces, moments, temperature,
see:ref:nbc <mdl.nbc>). If not specified with the `dof_ref`

parameter of the MDL nodes, the node-local base is
identical to the branch global base.

Note

Node related solutions are stored on the *b2000++* model database
in the branch-global coordinate, irrespective of nodal-local
transformations! See the laminate example case.

**Element Local Cartesian Base**

The element (local) coordinate system \(x_{el}, y_{el}, z_{el}\)
defines a base for computing material properties. *B2000++* makes the
following assumptions on how the element coordinate system is defined:

Beam elements: The element coordinate system is explicitly defined by the beam orientation.

Shell and continuum elements: The element coordinate system is implicitly defined by the first 3 element connectivity node identifiers \(n_1\), \(n_2\), and \(n_3\) defining the element, see also Generic Elements.

**Material Coordinate System**

Certain materials, like orthotropic or anisotropic materials, need to be formulated with respect to an own coordinate system called the material coordinate system \(x_{em}, y_{em}, z_{em}\). The mbase element property attribute defines the material orientation with respect to the current (branch) global orientation.

**Computational Coordinate System**

The computational coordinate system is the coordinate system in which
the computational degrees of freedoms of a node are expressed. The
term *computational* refers to the degrees of freedom which are used
in the final equation system and which express the relation between
displacements and applied forces, i.e. to the unknowns used in the
computational process.

On the node (DOF) level, the global computational coordinate system is usually identical with the node local coordinate systems. This is certainly the most convenient alternative for shell structures (with the exception of points located on a juncture line between several substructures). Also, for 3D structures the global computational and the node local systems usually coincide, or can be made to coincide with little difficulty.

*B2000++* requires two types of coordinate transformations, i.e
rotations (change of base) and translations. In what follows, we will
refer to two covariant coordinate systems, namely the master system
**m** described by base vectors \(\mathbf{m}_i\) and the slave
system **s** represented by base vectors \(\mathbf{s}_i\).

The change of base or rotation from \(\mathbf{s}\) to \(\mathbf{m}\) is defined by the scalar products of the base vectors \({m}_i\) and \({s}_i\):

A vector \(\mathbf{x}\) is rotated from \(\mathbf{s}\) to \(\mathbf{m}\) by the rotation matrix \(\mathbf{T}_{ij}\):

Note that the transformation is orthogonal. As an example, consider the transformation from the branch global coordinate system to the global coordinate system:

\(\mathbf{r}_{gg}\) designating the global (Cartesian) coordinate vector, \(\mathbf{r}_0\) the translation vector from the global coordinate system to the branch global coordinate system, and \(\mathbf{r}_{bg}\) the branch-global vector. \(\mathbf{T}\) is the rotation matrix described above.

All node related properties,i.e DOF’s, natural boundary conditions, or
essential boundary conditions are formulated with respect to the *node
local coordinate system* whose orientation is determined by a local
base relative to the branch coordinate system. Any node with a node
local coordinate system will have the transformation attribute specified by the `dof_ref`

option
of the node coordinate input specification and
contained in the third column of the NODA
dataset. By default, i.e if not specified otherwise, the node local
coordinate system is identical to the branch coordinate system. Any
node \(i\) which has a local coordinate system will have an index
\(k > 0\) in the NODA table. The index \(k\) then points to
the \(k^{th}\) row of the NLCS array
containing the rotation matrices.

## Loads in Nonlinear Analysis

All nonlinear analysis loads are by default scaled by the load
factor. This scaling can be altered by means of the `sfactor`

and
`sfunction`

directives in the case command.

**Example**: Gravitational forces with an acceleration of 3.5g are
applied to all elements:

```
nbc 1 type accelerations
accelerations 0 0 -9.806
allelements
end
case 1
analysis nonlinear
nbc 1 sfactor 3.5
end
```

Maintaining a constant loading throughout the analysis can be achieved
by means of the `sfunction 1`

option when specifying the `nbc`

set
identifier in the case command.

```
nbc 1 sfunction 1
```

A smooth load function can be defined in the `nbc`

set identifier in
the case command. The parameter `t`

is the built-in name for
the load factor (between 0 and 1 by default):

```
nbc 1 sfunction "0.5*(1-cos(PI*t))"
```

## Laminates

*B2000++* allows for the specification of laminated materials
(composite materials) for shell and - to a certain extent - continuum
elements. To work with laminated material a number of precautions have
to be taken, because:

Laminate layer orientation angles (see MDL laminate command) are always defined with respect to the current material reference coordinate system.

The default element material reference coordinate system is “branch global” or “global”.

A wrong element material reference coordinate system can therefore lead to wrong laminate orientations.

**Example**: A one-element test defines (1) an element with the
element-local coordinate system equal to the global coordinate
system (2) the element rotated by 45° around the global
z-axis. The element material is layered, with one layer with layer
orientation angle 0°. The element is statically determinate and
pulled along the right edge with a line load.

Results of model (2) are correct because of the correct mbase specification `mbase 12 0 0 0 0 0 0 0 0 0`

:

The material reference coordinate system is defined by angles of orientation (Euler angles) with respect to the element-local coordinate system (here: All Euler angles are 0, i.e that layer orientation is with respect to the element local system).

Note

Element stress and strain tensors are stored on the *B2000++* model
database with respect to the global coordinate system.