A thin quadratic laminated plate is loaded with a forced displacement along one of the plate edges. The linearised buckling load due to forced displacement of the edge is to be determined.
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 eltype Q9.S.MITC p1 0. 0. 0. p2 200. 0. 0. p3 200. 200. 0. p4 0. 200. 0. thickness 2 ne1 10 # N. of elements along edges 1 and 3 ne2 10 # N. of elements along edges 2 and 4 mid 2 # Material identifier end
Additional information passed to the patch is the material
identifier mid
describing the element material (see
below).
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. Layer 1 is the bottom layer and layer 16 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.
The material command defines the element material properties, grouping them by their material identifiers. In our case we have to define two materials, one for the laminate and one for the perlayer material:
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.e6 # tension, microstrain c 3000.e6 # compression, microstrain s +5000.e6 # shear, microstrain filter max_of_element end 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
Although this buckling load still within the elastic limit of the material, a simple failure criterion is used, based on the principals (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 (in this case by a factor of 16*27=432), 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) 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 the y and zdisplacement DOFs 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 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
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 1 value 1. dof FZ epatch 1 p9 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.
All parameters related to an analysis case description are contained in one command block, the case command block, which is referenced by the adir command:
case 1 analysis linearised_prebuckling nmodes 10 ebc 1 gradients 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
Results can also be viewed with the B2000++ browser. Enter the shell command below and click on "Solutions", "case 1", and "Buckling Loads":
$
b2browser buckling
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.
The buckling modes can be plotted with the baspl++ tool (B2000++ Pro). Ente the following shell command
$
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: