2. Buckling of Composite Plate

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 simply-supported on all four edges, one edge being free to move in the x-direction. The plate is loaded with a forced displacement loading condition along the free edge 2. For post-buckling analysis, a small out-of-plane load in the z-direction is applied at the center of the plate.

Geometry and loading of plate

Figure 7. Geometry and loading of 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 103 N/mm2, υ = 0.33

  • Laminate: E1 = 146.86 103 N/mm2, E2 = 9.65 103 N/mm2, E3 = 9.65e 103 N/mm2, υ12 = 0.3, υ23 = 0.023, υ13 = 0.023, G12 = 45.5 102 N/mm2, G23 = 45.5 102 N/mm2, G13 = 45.5 102 N/mm2

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 nine-node (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 counter-clockwise, 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

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 per-layer 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.e-6 # tension, micro-strain
    c -3000.e-6 # compression, micro-strain
    s +5000.e-6 # shear, micro-strain
    filter max_of_element

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

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 x-direction 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 z-displacement DOFs of all edge-nodes. The second line applies a x-displacement of -1 mm to the nodes of edge 2, while the third line locks the x-displacement 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

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 out-of-plane force at the center node, is defined with the nbc command:

nbc 1
  value 1. dof FZ epatch 1 p9

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

  case 1

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.1797e-02 +6.6287e+03
   2 +9.2413e-02 +1.4656e+04
   3 +1.1776e-01 +1.8675e+04
   4 +1.5845e-01 +2.5130e+04
   5 +1.8564e-01 +2.9441e+04
   6 +2.3845e-01 +3.7817e+04
   7 +2.4416e-01 +3.8723e+04
   8 +2.8847e-01 +4.5750e+04
   9 +3.0698e-01 +4.8685e+04
  10 +3.3667e-01 +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 re-running 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 post-buckling path.

The buckling modes can be plotted with the baspl++ tool (B2000++ Pro). Ente the following shell command

$ baspl++ view-buckling.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: