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
and7
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 incrementaliterative algorithm, also called predictorcorrector 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 firstorder variation.
Thus we have:
To impose the constraint, define Lagrange multipliers \(\Lambda\) to form the Lagrangian:
The firstorder derivatives are:
where \(c\left( v,\Lambda,\lambda \right)\) is the increment constraint condition that defines the increment control strategy.
The secondorder derivatives are:
Nonlinear Dynamic Solver Numerical Algorithm
The dynamic nonlinear problem is:
The general form of a LMS (linear multistep) scheme for a normalized firstorder ODE (ordinary differential equation) of the form:
is:
where
and
Adapted to a secondorder ODE (ordinarydifferential equation) on u where:
we obtain:
where
At each time integration step, the problem to solve is:
The firstderivative of the potential energy of the unconstrained problem is:
To impose the constraint, add Lagrange multipliers to form the Lagrangian:
The firstorder derivatives are:
The secondorder derivatives are:
Linearised PreBuckling 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:
FrequencyDependent Vibration Solver Numerical Algorithm
The frequencydependent 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 FrequencyDependent Undamped FreeVibration 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.0e4). 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 freevibration
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 wellknown problem with shell elements is how intersections or junctions (such as the connecting edge between panels and stiffeners) are treated. Pure 5 degreesoffreedom (DOF) elements are not able to take into account the drill moment and the rotation in the element outofplane direction (local zdirection). This violation of the compatibility requirement results in numerical problems that often prevent a successful simulation. Elements with an additional degreeoffreedom 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 drillstabilized. This is done such that junction nodes and nodes where boundary conditions or linear multipoint 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++ postprocessor 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 degreeoffreedom 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
nonsymmetric 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 metainformation to the database. The B2000++ solvers produce 3 types of solution fields:
DOF fields: DOF fields are 2dimensional arrays containing solution data define at the nodes 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 are 2dimensional arrays containing solution data or derived data defined at the elements of the mesh.
Sampling point fields are elementwise derived quantities of the DOF fields. The values are stored elementwise 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.
All solution fields are identified by the following parameters identifying the field on the database:
gname.branch.cycle.subcycle.case.subcasemode
Parameters for identifying fields with tools, such as Simples or baspl++ are equal to those listed below.
Parameter 
Description 


Generic name, such as 

Branch or mesh identifier, a positive integer. 

The computational cycle identifier, a positive integer or 0. 

The computational subcycle identifier, a positive integer or 0. 

The analysis case (‘load’ case) identifier,a positive integer. 

The analysis subcase identifier, a positive integer. Subcases identify linear analysis case with common constraints or eigenmodes originating from eigenvalue analysis. 
Deformation Analysis Solution Fields
For static and transient analysis the solvers of B2000++ produce the following solution datasets on the model database:
Name 
Description 


Displacement and rotation fields. Fields usually defined at the mesh nodes. processor. Generated by the solver. 

Velocity fields (transient analysis). Fields usually defined at the mesh nodes. Generated by the solver. 

Acceleration fields (transient analysis). Fields usually defined at the mesh nodes. Generated by the solver. 

Forces and optional moments (usually) applied to the mesh nodes. Generated by the solver. 

‘Reaction’ forces and moments computed by the relevant solver (usually) at the mesh nodes. 

Vibration modes (normalized displacement and rotation fields).The descriptors of these sets contain the solution attributes, i.e circular frequency, frequency, modal stiffness, modal mass. Generated by the solver. 

Buckling modes (normalized displacement and rotation fields). The descriptors of these sets contain the solution attributes, i.e eigenvalue, critical load. Generated by the solver. 
Name 
Description 


Sampling point fields containing strain tensors with 6 components at integration points. The reference system is global. Generated by the gradient solver. 

Sampling point fields containing stress tensors with 6 components at integration points. The reference system is global. Generated by the gradient solver. 

Sampling point fields containing strain tensors with 1 component at rod and cable element sections. The reference system is elementlocal. Generated by the gradient solver. 

Sampling point fields containing stress tensors with 1 component at rod and cable element sections. The reference system is elementlocal. Generated by the gradient solver. 

Sampling point fields containing 3 force components at rod element sections. The forces are computed with respect to the deformed system (nonlinear analysis). The reference system is global. Generated by the gradient solver. 

Sampling point fields containing 3 force components at beam element sections. The reference system is global. The forces are computed with respect to the deformed system (nonlinear analysis). Note that stress calculations for beams is not integrated in the B2000++ solver. The Simples bc module computes section stresses. Generated by the gradient solver. 

Sampling point fields containing 1 force component at rod element sections in the rod direction (local xaxis). Generated by the gradient solver 

Sampling point fields containing 3 moment components at beam element sections. The reference system is global. Note that stress calculations for beams is not integrated in the B2000++ solver. The Simples bcs module computes section stresses. Generated by the gradient solver. 
Heat Flow Analysis Solution Fields
The heat flow analysis solver of B2000++ produces the following solution datasets on the model database:
Name 
Description 


Concentrated heat values acting (usually) at the mesh nodes. Generated by solver. 

Concentrated heat values computed by the relevant solver (usually) at the mesh nodes. 

Temperatures computed by the relevant solver (usually) at the mesh nodes. Generated by the gradient solver 

Heat flux at element sampling points. Generated by the gradient solver 

Heat capacity at element sampling points. Generated by the gradient solver 
Field Attributes
The FIELDS dataset contains dictionaries describing attributes of all B2000++ solution fields. They help postprocessors and other external programs to decide what to do with a particular field.
Key 
Value 

GNAME 
FORC 
ADDITIVE 
YES or NO 
DISCRETE 
YES or NO 
INDEXED 
YES or NO 
SYSTEM 
BRANCH 
TYPE 
NODE or ELEMENT 
Coordinate Systems and Transformations
B2000++ is substructure (branch) oriented. A B2000++ mesh can consist of an arbitrary assemblage of meshes collected in branches connected by nodes. However, these meshes or branches themselves may not contain other branches.
B2000++ defines several Cartesian coordinate systems for describing the orientation of nodes, elements, materials, DOFs, stresses, etc.:
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 nodelocal base is
identical to the branch global base.
Note
Node related solutions are stored on the B2000++ model database in the branchglobal coordinate, irrespective of nodallocal 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 branchglobal 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 builtin name for
the load factor (between 0 and 1 by default):
nbc 1 sfunction "0.5*(1cos(PI*t))"
Laminates
B2000++ supports 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 oneelement test defines (1) an element with the elementlocal coordinate system equal to the global coordinate system (2) the element rotated by 45° around the global zaxis. 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 elementlocal 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.
Mass and nonstructural mass
B2000++ defines the non structural mass of rod, beam, ans shell elements
in the property block for beam elements (Bx.S.RS). Reason: The Bx beam elements are the only B2000++ elements the define the element properties with the property block  similar to BDF, referring to pid (Nastran BDF defines the non structural mass in the element property cards, not in the element cards). Example definition nonstructural mass:
property 11 type beam_section mid 1 shape rectangle d1 0.01 d2 0.01 non_structural_mass 0.27 end
in the element blocks for all other elements concerned. Reason: property blocks are not understood by these elements. Example definition nonstructural mass:
elements ... eltype Q4.S.MITC non_structural_mass 0.27 ... end