10. Buckling and Post-Buckling Analysis of Plate

A thin quadratic laminated plate is loaded with a forced displacement along one of the plate edges. The linearised buckling behaviour and the post-buckling 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++ view-buckling.mdl, view-postbuckling.py

10.1. Specifications

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

10.2. The Finite Element Model

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
  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 per-layer 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 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 DOFs the y- and z-displacements 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
  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 out-of-plane 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 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.

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.

10.3. Buckling Analysis

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

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.

Using the following Python script view-buckling.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++ 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:

10.4. Post-buckling Analysis

While the non-linear 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 pre-buckling 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 post-buckling 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 load-end 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 non-linear analysis options.

The post-buckling 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 view-postbuckling.py is written for the present example and extracts several XY curves:

$ baspl++ view-postbuckling.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 view-postbuckling.py works also with several databases. This allows to compare different discretizations, materials, boundary conditions, solution parameters, etc. In the file postbuckling-fast.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++ postbuckling-fast.mdl

and invoking the script

$ baspl++ view-postbuckling.py postbuckling.b2m postbuckling-fast.b2m

shows that, if the step size is not restricted, the unstable path is followed and no buckling is detected.