A thin quadratic laminated plate is loaded with a forced displacement along one of the plate edges. The linearised buckling behaviour and the postbuckling behaviour have to be studied. Some of the basic capabilities of B2000++, i.e solving nonlinear structural problems are demonstrated. The following files are found in the test directory:
MDL files  plate.mdl , buckling.mdl ,
postbuckling.mdl 
Python scripts for baspl++  viewbuckling.mdl ,
viewpostbuckling.py 
The plate to be analyzed is quadratic 200 mm by 200 mm with a thickness of 2 mm. It is simplysupported on all four edges, one edge being free to move in the xdirection. The plate is loaded with a forced displacement loading condition along the free edge 2. For postbuckling analysis, a small outofplane load in the zdirection is applied at the center of the plate.
Two materials are studied, an isotropic one (aluminum) and a layered composite material, 16 layers of 0.125mm thickness, with the symmetric stacking sequence (+45°,0°,45°,0°,+45°,0°,45°,+90°)_{s}. The material data are listed below:

Aluminum: E = 69 10^{3} N/mm^{2}, υ = 0.33

Laminate: E_{1} = 146.86 10^{3} N/mm^{2}, E_{2} = 9.65 10^{3} N/mm^{2}, E_{3} = 9.65e 10^{3} N/mm^{2}, υ_{12} = 0.3, υ_{23} = 0.023, υ_{13} = 0.023, G_{12} = 45.5 10^{2} N/mm^{2}, G_{23} = 45.5 10^{2} N/mm^{2}, G_{13} = 45.5 10^{2} N/mm^{2}
Now we explain in detail how to get to the solution. We
start by defining the Finite Element model by means of the
B2000++ model description language
(MDL). MDL contains some simple mesh generation capabilities by
means of epatch (element patches) command. A
single patch is used, the patch geometry type
plate
is selected, and a regular mesh of 10 by 10
elements is generated, selecting the Q9.S.MITC
element, i.e the ninenode (Q9) shell (S) element based on the
MITC formulation. The plate geometry is defined by the four
corner points p1, p2, p3, and p4 specified counterclockwise, as
well as the thickness (note that the thickness is ignored in case
of laminates, since the laminate defines the thickness). The
corresponding MDL block is listed below:
epatch 1 geometry plate p1 0. 0. 0. p2 200. 0. 0. p3 200. 200. 0. p4 0. 200. 0. thickness 2 eltype Q9.S.MITC mid 2 # Material identifier ne1 10 # N. of elements along edges 1 and 3 ne2 10 # N. of elements along edges 2 and 4 end
Additional information passed to the patch is the material
identifier mid
describing the element material
(see below).
The material command defines the element material properties. In our case we have to define three materials, on for the perlayer material, one for the laminate, and one for the isotropic material:
material 1 type orthotropic e1 146.86e+03 e2 9.65e+03 e3 9.65e+03 nu12 0.3 nu23 0.023 nu13 0.023 g12 45.5e+02 g23 45.5e+02 g13 45.5e+02 end material 2 type laminate 0.125 45 1 0.125 0 1 0.125 45 1 0.125 0 1 0.125 45 1 0.125 0 1 0.125 45 1 0.125 90 1 0.125 90 1 0.125 45 1 0.125 0 1 0.125 45 1 0.125 0 1 0.125 45 1 0.125 0 1 0.125 45 1 end material 3 type isotropic e 69e+03 nu 0.33 end
The laminate is defined by listing all layers one after the other, each layer specifying the layer thickness, the orientation of the layer with respect to the selected material identifier, and the material identifier associated to the layer. The first layer is the bottom layer and the last layer is the top layer. Note that, 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.
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) for creating the required displacement constraints along edge 2 in the xdirection and (2) locking the structure along the edges.
Essential boundary conditions (ebc) along patch edges, i.e displacement constraints, can be added in the patch. The patch defines for each edge a list of nodes which can be identified with e.g. epatch 1 e2, in this context meaning all nodes of edge 2 of epatch 1.
Since the edges are free to rotate we only lock the displacement DOFs, i.e DOFs UX, UY, or UZ. The first line locks DOFs the y and zdisplacements of all edgenodes. The second line applies a xdisplacement of 1 mm to the nodes of edge 2, while the third line locks the xdisplacement DOF of all nodes belonging to edge 4.
ebc set 1 value 0. dof [ UY UZ] epatch 1 e1 epatch 1 e2 epatch 1 e3 epatch 1 e4 value 1. dof [UX ] epatch 1 e2 value 0. dof [UX ] epatch 1 e4 end end
The ebc set will be added to the boundary condition sets of the case to be treated by the analysis, as will be described later.
The natural boundary condition (nbc) of our model, i.e the outofplane force at the center node, is defined with the nbc command:
nbc set 1 value 1. dof FZ epatch 1 p9 end end
A concentrated force of 1 N is assigned to the shell patch center node (this node is identified with epatch 1 p9).
Note that, if nodelocal 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 nodelocal coordinate systems, hence, all boundary conditions are expressed in the global coordinate system.
The elements generated by the patch then refer either to the
laminate material (mid 2
) or the isotropic
material (mid 3
). In the following, only the
laminate material is being used.
All this is combined in the file
plate.mdl
which is included by the files
buckling.mdl
and
postbuckling.mdl
(see below).
This terminates the description for the Finite Element model. What remains now is to add all the components defined so far in a description of the analysis case and to the tell the solver which analysis cases will be solved.
All parameters related to an analysis case description are
contained in one command block, the case command block. The file
buckling.mdl
contains these MDL
commands and the analysis directives as well:
include 'plate.mdl' case 1 analysis linearised_prebuckling nmodes 10 ebc 1 end adir case 1 end
This file includes the previouly described file
plate.mdl
. In the analysis case description,
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) are activated.
The buckling analysis is carried out by entering the shell command
$
b2000++ buckling.mdl
After the analysis has completed, the following command prints the buckling factors and the buckling loads:
$
b2print_modes prec=5 buckling.b2m
Linear reaction forces +1.5859e+05 Number of computed buckling modes for case 1: 10 Mode Load factor Critical force 1 +4.1797e02 +6.6287e+03 2 +9.2413e02 +1.4656e+04 3 +1.1776e01 +1.8675e+04 4 +1.5845e01 +2.5130e+04 5 +1.8564e01 +2.9441e+04 6 +2.3845e01 +3.7817e+04 7 +2.4416e01 +3.8723e+04 8 +2.8847e01 +4.5750e+04 9 +3.0698e01 +4.8685e+04 10 +3.3667e01 +5.3393e+04
A critical load of 6.63 kN is predicted. Changing the values
for ne1 and ne2 to 20 in the
file plate.mdl
, and rerunning the analysis,
one obtains a critical load of 6.59_kN, a difference of
0.6%. Hence, the 10x10 mesh is sufficiently fine to predict the
first buckling and the initial postbuckling path.
Using the following Python script
viewbuckling.py
m = Model(sys.argv[1]) p = NPart(m) p.contour.colmap.texture.mode = 'step' p.elements.extract = True o = Modes() o.part = p
and entering the following command on the shell
$
baspl++ viewbuckling.py buckling.b2m
displays the plate with the first buckling mode. Entering the command
>>>
o.field.ch.mode = 2
on the Python prompt selects the second buckling mode. The following image shows the shape's of the buckling modes, only the first buckling mode is of interest:
While the nonlinear analysis case description refers to the
same discretization, materials, etc. as the buckling analysis, the
solution strategy is different, which is shown in the file
postbuckling.mdl
:
include 'plate.mdl' case 1 analysis nonlinear ebc 1 sfactor 0.1 # up to 0.1mm edge shortening nbc 1 sfunction '0.1' # constant force, not scaled by load factor step_size_init 0.05 # to get some points on prebuckling path step_size_max 0.05 # to get some points gradients 1 # compute gradients for all increments end adir case 1 end
The essential boundary condition set ebc 1 must be scaled; from the linearised buckling analysis it is known that first buckling occurs at 0.042 mm edge shortening. In this example, the postbuckling analysis should go up to 0.1 mm edge shortening; this is achieved with the sfactor command. Note that sfactor 0.1 is identical to sfunction '0.1*t', where t is the current load factor or time.
The concentrated force constitutes an artificial imperfection to leave the unstable path and jump to the stable path. It is not scaled by the load factor but fixed to 0.1 N (by using with the sfunction '0.1' command).
A rather delicate choice is the step size. The step sizes to be selected depend on what is desired: To find, for this simple example, the stable path and to obtain a loadend shortening curve with enough points, the maximum step size (here: 0.05) has to be set, producing a minimum of 1/0.05 = 20 increments. The initial step size is also crucial for a successful analysis: A good choice is a fraction of the buckling load. See here for more information on static nonlinear analysis options.
The postbuckling analysis is carried out by entering the shell command
$
b2000++ postbuckling.mdl
To monitor the progress of the analysis, one can enter instead:
$
b2000++ l 'info of solver in cerr' postbuckling.mdl
The Python script viewpostbuckling.py
is written for the present example and extracts several
XY curves:
$
baspl++ viewpostbuckling.py postbuckling.b2m
The strain curve shows best where buckling occurs, at 6.50 kN and 0.041 mm, which is very close to the value that was predicted by the buckling analysis.
The script viewpostbuckling.py
works
also with several databases. This allows to compare different
discretizations, materials, boundary conditions, solution
parameters, etc. In the file
postbucklingfast.mdl
the step size is not
restricted:
include 'plate.mdl' cases case 1 analysis nonlinear ebc 1 sfactor 0.1 # up to 0.1mm edge shortening nbc 1 sfunction '0.1' # constant force, not scaled by load factor gradients 1 # compute gradients for all increments end end adir case 1 end
Analysis of this file
$
b2000++ postbucklingfast.mdl
and invoking the script
$
baspl++ viewpostbuckling.py postbuckling.b2m postbucklingfast.b2m
shows that, if the step size is not restricted, the unstable path is followed and no buckling is detected.