# 3. Buckling Analysis

Folder: $B2EXAMPLES/buckling/cylinder_cutout A cylinder with two opposite square cutouts is clamped at one end and compressed at the other end in the axial direction by an imposed displacement of 0.00425. This problem has been proposed by Almroth and Holmes [almroth1972] and has later been taken up by Knight and Rankin [knight2006]. The goal of this example is to demonstrate the collapse analysis and the post-processing capabilities of B2000++. Note that the B2000++ Verification Manual treats this case, too, but with Q4 and Q9 elements. Similarly to Knight and Rankin, the problem is modeled with an eight of a cylinder (half in the longitudinal direction and a fourth in the circumferential direction): To this end, 2 cylindrical epatch patches are defined. The figure below displays the mesh with Q4.S.MITC elements. Boundary conditions are symmetry along both longitudinal edges (edges defined by blue points P1-P4 and green points P2-P3) and along one circumferential edge (edges defined by blue points P3-P4), the other circumferential edge being clamped. The test is run with Q4.S.MITC.E4 shell elements, comparing the cutout tip displacement-load function (point A) against the reported solution by Knight and Rankin. The analysis control block contains a standard set of options, the non-linear analysis method being the default continuation method): case 1 analysis nonlinear ebc 1 ebc 2 step_size_init 0.05 step_size_max 0.05 gradients 1 rcfo_restrict nodeset "EndShortening" end  The maximum step size limitation is imposed in order to obtain a load curve with some points. Note the rcfo_restrict option: It contains the node set(s) containing all nodes involved in the imposed displacement, thus rendering post-processing very easy (see stroke-load diagram below). The y-axis of both diagrams below are scaled to the reaction force $$P/P_{0}$$, with $$P_{0} = {{1.2}\pi Et^{2}}$$, where $$E=1 \times 10^{7}$$ and $$t=0.014$$. The Simples script run.py (1) runs the analysis for selected mesh densities, extracting the graphs for each analysis and (2) collects the graphs and produces the stroke versus load disp_load.png and stroke_load.png out of plane versus load plots. # Python script for executing a series of B2000++ runs with increasing # mesh density for different element types, followed by a synthesis # graph. import os import simples as si # Comparison values obtained with STAGS stags_dx = [0.00891496, 0.0137947, 0.0204575, 0.0290596, 0.0404457, 0.0551476, 0.0731026, 0.0914018, 0.0975327, 0.0991906, 0.0989717, 0.102444, 0.105009, 0.107511, 0.109607, 0.112328] stags_load = [0.086217, 0.104633, 0.119883, 0.137595, 0.162698, 0.205865, 0.277067, 0.365396, 0.373607, 0.375367, 0.379707, 0.386628, 0.394956, 0.395191, 0.38393, 0.376657] os.system("rm -rf *.b2m *.modeler") # Compute P0 = Analytical buckling load of cylinder E = 1e7 t = 0.014 P0 = 0.25*1.2*3.14159*E*t**2 # Run analysis for different mesh densities for one element type eltype = 'Q4.S.MITC.E4' radial_disp_reaction_graphs = [] stroke_reaction_graphs = [] meshes = ((18,32), (36,64)) for ne in meshes: mesd = '%s %dx%d' % (eltype, ne[0], ne[1]) print('Run for mesh density %s' % mesd) # Run B2000++ for current parameters if si.util.solve('cutout.mdl', params={'ELTYPE': eltype, 'NE1': ne[0], 'NE2': ne[1]}): # Open db and get node id of point A. model = si.Model('cutout.b2m', 'or') # Get internal node number of epatch 1 point P3 == point A. point = model.get_node_set('EPATCH-1-P')[2] # Get all analysis steps values. Scale because TIME goes from # 0 to 1, but displacement constrains from 0 to 0.1 steps = model.get_steps(scale = 0.1) # Get the displacements of node *point*. Extract the # node-local component 'radial', which is dof 1. dy = model.get_dof_value_of_steps('displacements', node = point, dof = 1, system = 'nl') # Get sum of all reaction forces due to displacement loading and scale. ry = model.get_reaction_force_of_steps(scale = 1./P0) model.close() # Create Graph for current model and add to list radial_disp_reaction_graphs.append((mesd, (dy, ry))) stroke_reaction_graphs.append((mesd, (steps, ry))) # Plot all unctions colors = ('red', 'green', 'blue', 'cyan', 'magenta', 'yellow') markerw = 0.5 markers = 3. g = si.Graph2() for i, a in enumerate(radial_disp_reaction_graphs): xy = a[1] line = g.add_line(xy[0], xy[1]) line.marker= 'o' line.markeredgewidth = markerw line.markersize = markers line.markerfacecolor = None line.color = colors[i] line.label = a[0] line = g.add_line(stags_dx, stags_load) line.color = 'black' line.label = 'STAGS' g.axes_label_x = 'Radial displacement at A' g.axes_label_y = 'Load P/P0' g.add_legend((.2, .7)) g.save('disp_load.png') g = si.Graph2() for i, a in enumerate(stroke_reaction_graphs): xy = a[1] line = g.add_line(xy[0], xy[1]) line.marker= 'o' line.markeredgewidth = markerw line.markersize = markers line.markerfacecolor = None line.color = colors[i] line.label = a[0] g.axes_label_x = 'Stroke' g.axes_label_y = 'Load P/P0' g.add_legend((.2, .7)) g.save('stroke_load.png')  ## 3.2. Buckling of Three Stringer Composite Panel Folder: $B2EXAMPLES/buckling/panel3

The three string composite panel is a test case for demonstrating the B2000++ capabilities of dealing with lightweight composite structures modeled with shell or solid elements. A cylindrical panel with 3 stiffeners is subjected to axial compression load. The lowest buckling loads are to be determined.

The boundary conditions are displayed in the figure below: Both curved edges are clamped and both straight edges are free. One of the curved edges is pushed towards the panel, causing it to buckle. All boundary conditions are essential boundary conditions (ebc command in the MDL input files) Thus, the resulting reaction force due to the edge displacement in axial direction is equal to the buckling load when multiplied by the computed eigenvalue.

The shell model consists of Q9 MITC shell elements. The cylindrical shell and the shell plus the stringer foot are modeled with one element through the thickness, as is the stringer. The laminate orientations and laminate placement are visualized with baspl++:

Buckling modes without prestress of the shell model are computed for load case 1, requiring the first 10 modes:

case 1
ebc       1
analysis  linearised_prebuckling
nmodes    10
end


The buckling analysis is launched by the shell command

prompt> b2000++ --adir 'cases=1' -l "info of solver in cerr" shell.mdl


Linear buckling analysis results can be extracted and printed with b2print_buckling_loads:

b2print_buckling_loads shell

Linear reaction forces +2.9477108e+05

Number of computed buckling modes for case 1: 10
1 +1.7227278e-01 +5.0781032e+04
2 +1.7237261e-01 +5.0810459e+04
3 +1.7464288e-01 +5.1479669e+04
4 +1.7469178e-01 +5.1494083e+04
5 +1.7908786e-01 +5.2789922e+04
6 +1.7965825e-01 +5.2958056e+04
7 +1.8222893e-01 +5.3715819e+04
8 +1.8233669e-01 +5.3747582e+04
9 +1.8770931e-01 +5.5331277e+04
10 +1.8883024e-01 +5.5661694e+04


The first buckling mode shape plot is obtained with the baspl++ script view_buckling_modes.py.

Folder: $B2EXAMPLES/buckling/plate A simply-supported thin quadratic laminated plate of 500 mm by 500 mm is loaded with a forced displacement along one of the plate edges. The linearized buckling load as well as the post-buckling response due to the forced displacement of the edge is to be determined. The plate is simply supported on all four edges, one edge being free to move in the x-direction to allow for the forced displacement. The plate is loaded with a forced displacement loading condition along the free edge 2. The plate is made of a composite material, 16 layers of 0.125mm thickness, with the symmetric stacking sequence (in degrees) $$(+45,0,-45,0,+4,0,-45,+90)_s$$. The material data are listed below: $$E_1 = 146.86\ 10^{3} N/mm^2$$ $$E_2 = E_3 = 9.65\ 10\ 10^{3} N/mm^2$$ $$\nu_232 = \nu_13 = 0.023$$ $$G_12 = G_23 = G_13 = 45.5\ 10^{2} N/mm^2$$ FE Model Now we explain in detail how to get to the solution. We start by defining the Finite Element model with the B2000++ model description language (MDL). MDL contains simple mesh generation capabilities with element patches. A single epatch defines the regular plate mesh and links to the material. (et="Q4.S.MITC.E4") (ne=40) # Mesh definition epatch 1 geometry plate # Geometry type p1 0. 0. 0. # 4 plate vertices p2 500. 0. 0. p3 500. 500. 0. p4 0. 500. 0. eltype (et) # Element type ne1 (ne) # N. of elements in first direction ne2 (ne) # N. of elements in second direction mid 2 # Link to material identifier end  Note the parametrization of the element type and the number of elements. Here we select Q4 four node quadrilateral MITC elements. The plate geometry is uniquely defined by the four corner points p1, p2, p3, and p4 specified counter-clockwise, defining a flat plate. Note that the plate thickness is not needed if the material is a laminate material. A laminate is defined with the MDL material block, by listing all layers one after the other, each layer specifying the thickness, the layer orientation, and the material identifier. # Composite stacking sequence definition material 2 type laminate 0.125 45 1 0.125 0 1 ... end  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 composite material constants are defined with a MDL material block, specifying the orthotropic material type. 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 end end  A simple failure criterion is used, based on the maximum of the principal 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, 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) creating the required displacement constraints along edge 2 in the x-direction and (2) properly locking the structure along the edges. Essential boundary conditions ebc along patch edges, i.e displacement constraints, can be added by referencing the proper epatch element (here: edge e2). MDL automatically generates the mesh node and element lists. They are stored on the model database and can be referenced by ebc and nbc boundary condition definitions. The first ebc block 1 locks the structure. Since the edges are free to rotate we only lock the displacement DOFs, i.e DOFs UX, UY, or UZ:: ebc 1 # Zero-value constraints value 0. dof [ UY UZ] epatch 1 e1 epatch 1 e2 epatch 1 e3 epatch 1 e4 value 0. dof [UX ] epatch 1 e4 end  The first line of the ebc block locks the y- and z-displacement DOFs of all edge nodes. The second line locks the x-displacement DOFs of the edge opposite to the moving edge (e2). The second ebc block 2 applies a value for the x-displacement of all nodes of the edge e2: ebc 2 # Forced displacement value -1. dof UX epatch 1 e2 end  The links to the ebc sets will be added to the case block, as will be described later. For the particular configuration, i.e. a flat plate, an ‘imperfection’ has to be defined to trigger the post-buckling behavior. This can be achieved by defining a small out-ouf-plane force that ‘pushes’ the plate in the desired post-buckling direction (in the present model only one configuration is possible). The out-ouf-plane concentrated force of 1 N is assigned to the shell patch center node (identified with the epatch node p9) with a nbc block: # Central out-of-plane force to create imperfection nbc set 1 value 1. dof FZ epatch 1 p9 end end  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 a case block. In out example we will carry out 2 types of analysis, a buckling analysis and a post-buckling analysis. Create model database b2ip++ plate.mdl  Buckling Analysis The case block for the buckling analysis defines the type of analysis, includes the links to the boundary conditions, and, for the buckling analysis type, defines the number of buckling modes to be computed. case 1 analysis linearised_prebuckling nmodes 10 ebc 1 ebc 2 gradients 1 end adir # needs at least one case for -adir cases= to work case 1 end  In the analysis case description above, 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 and ebc 2) are activated. The buckling analysis (case 1) is carried out by entering the shell command prompt> b2000++ -adir cases=1 plate.b2m  After the analysis has completed, the b2print_buckling_loads prints the critical loads: b2print_buckling_loads buckling Eigenvalues and buckling loads, case 1 Contributing EBC identifiers: 2 Total reaction force F=(-158593, 2.50111e-12, -3.66168e-15) amplitude Famp=158593 Mode Lambda Crit. load 1 0.01718 2725 * 2 0.03948 6262 * 3 0.05028 7973 * 4 0.06912 1.096e+04 * 5 0.08614 1.366e+04 * 6 0.1099 1.743e+04 * 7 0.1128 1.789e+04 * 8 0.1335 2.117e+04 * 9 0.1628 2.582e+04 * 10 0.1632 2.588e+04 * Lambda (eigenvalue): if > 1 no buckling, < 1 buckling, marked by * Critical load: F * Lambda (if F != 0) Reaction force: Sum of reactions due to EBC and NBC  Results can also be viewed with the b2browser. Enter the shell command below and click on “Solutions”, “case 1”, and “Buckling Loads”: prompt> b2browser buckling  The buckling modes can be plotted with the baspl++ viewer. Enter the following shell command which displays the plate with the first buckling mode: prompt> baspl++ plot_buckling_mode.py  The following image shows the shapes of some buckling modes: Post-Buckling Analysis The case block for the post-buckling analysis case 2 defines the type of analysis, includes the links to the boundary conditions, and, for the buckling analysis type, defines the number of buckling modes to be computed. case 2 analysis nonlinear ebc 1 ebc 2 sfactor 0.1 # up to 0.1mm edge shortening nbc 1 sfunction 1. # constant force, not scaled by load factor step_size_init 0.01 # to get smoooth pre-buckling path step_size_max 0.01 # to get smoooth post-buckling path gradients 1 # compute gradients for all increments rcfo_restrict epatch 1 e2 # Restrict computation of reaction forces # to edge e2 nodes end  The following comments pertain to the the post-buckling analysis case description above: • The ebc 2 link includes the sfactor 0.1, which means that the load is incremented up to 0.1 times the ebc 2 values. • The nbc 1 link includes sfunction 1.0, which means that the load of nbc 1 is applied constant during all increments. • Setting step_size_init and step_size_max both to 0.01 , i.e. keeping the increments constant, will produce a smooth post-buckling path for the graphs, at the expense of higher computation time.2 • The rcfo_restrict option tells any post-processing tool, such as Simples, to restrict the computation of reaction forces to the nodes of the list (here: the epatch edge e2 nodes). Post-buckling analysis is carried out by entering the shell command prompt> b2000++ -adir cases=2 plate.b2m INFO:solver.static_nonlinear:18:22:45.812: Start static nonlinear solver for case 1 INFO:solver:18:22:45.813: Resolving linear constraints: 102 equations. INFO:solver:18:22:45.813: 0 dependent constraints were eliminated. INFO:solver:18:22:45.813: Reduction matrix (903 x 1005) has 903 nonzeros. INFO:solver.static_nonlinear:18:22:45.813: Start increment 1 of stage 1 for time 0 and time step 0.01. INFO:solver.static_nonlinear:18:22:45.885: Start increment 2 of stage 1 for time 0.01 and time step 0.01. ... INFO:solver.static_nonlinear:18:22:48.304: Start increment 99 of stage 1 for time 0.98 and time step 0.01. INFO:solver.static_nonlinear:18:22:48.331: Start increment 100 of stage 1 for time 0.99 and time step 0.01. INFO:solver.static_nonlinear:18:22:48.363: End static nonlinear analysis for case 1 with success.  After the analysis has completed, the simples script plot_postbuckling_graphs.py prompt> python plot_post-buckling_graphs.py  plots the graphs below: ## 3.4. Restart of Curved Panel Folder: $B2EXAMPLES/buckling/restart

The test case computes the post-buckling behavior of a slightly curved panel, thereby testing the restart capabilities of B2000++.

The panel has a length L=92mm along the straight edges, a radius R=586mm with an angle α=9.0 degrees along the curved edges and a thickness t=1mm. The material is linear isotropic, with modulus of elasticity E=70000 Pa and Poisson number 0.3. The edge loads act on edge I in the positive z-direction. Displacements in radial direction on edge I are suppressed and edge IV is simply supported. Edge II and III have symmetric boundary conditions. The mesh density is 6 by 6 Q4 elements.

Due to symmetry, only one quarter of the panel is considered. The panel is modeled with Q4.S.MITC.E4 shell elements.

The restart information is added to the database with the restart function (see file restart.py):

def restart(db, case=1, step=None, sfactor=2.0):
"""Modifies the database for an imminent restart from implicitly
assumed last step (step=None) or from specified step (cycle).

:param db db: pymemcom database object.

:param int case: Case identifier.

param int step: Step (cycle) from which to start analysis
again. Default is last step found.

:param foar sfactor: New sfactor to e reached at last step.

"""

stage = 1 # first stage of case 1
if not step:
laststep = db[f'SOLUTION-STAGE.0.0.0.{case}.{stage}'][:]['LAST_STEP'][0]
else:
laststep = step

# print(f'Restart from last step={laststep}, time={time}')
db[f'CASE.{case}'][0]['SFACTOR']=sfactor
db[f'CASE.{case}'].desc['DOF_INIT']=f'DISP.*.{laststep}.0.{case}'
db[f'CASE.{case}'].desc['STEP_INIT_NUM']=laststep
# time = db[f'SOLUTION-STEP.0.{laststep}.0.1'][:]['TIME'][0]
# db[f'CASE.{case}'].desc['TIME']=time
db.close()


## 3.5. Rod Snap Problem

Folder: \$B2EXAMPLES/buckling/snap

This simple problem consisting of two rods illustrates the snap phenomenon and its solution with B2000++. Note that this model is purely theoretical, because the material is assumed linear elastic!

The FE model is shown below, the rods being inclined by 10 degrees. Material is linear isotropic with a modulus of elasticity roughly corresponding to neoprene or similar. The geometric and material values can be consulted in file rod.mdl.

The model is constrained at nodes 1 and 3, node 2 being free to move in the x- and y-direction. A load (nbc) $$F_y=-200$$ is applied at node 2.

When trying to solve the problem for the natural boundary condition case, i.e the case with a force $$F_y$$ applied in the y-direction of node 2, the Newton-type solution procedure will fail, because the procedure will not overcome the situation where $$U_y=H$$. Arch-length methods are then indicated (see input files).

The applied load $$F_y$$ (figure above) of load case 1 is increased up to a factor 10. Since we do not increment from 0 to 1 we need to define at least one stage.

Two models are compared:

• A model with 1 stage with a SFACTOR of 10 (1stage.mdl). The case 1 block then adds a single stage (i.e case) 11, defining the constraint (ebc) set, the load (nbc) set, and the single stage 11:

case 1
analysis               nonlinear
increment_control_type hyperplane
stage                  11 sfactor 10.0

• A model with 2 stages, each with a SFACTOR of 5, leading to a final SFACTOR of 10 (2stages.mdl). Observe that, with the selected step increment parameters, the solution jumps to secondary branch! The case 1 block then adds 2 stages (i.e case) 11 and 12, defining the constraint (ebc) set, the load (nbc) set, and the stages 11 and 12:
case 1