Nonlinear Axial Compression (Buckling) of a Cylinder not Converging

[fethio]

thanks @Calc_em . So, S4R expands to a single layer of C3D8R, which suggests that it will be advisable to switch to a volume model and create more than one layer in the thickness direction, right? This is probably the reason why others have repeatedly advised the use of continuum solid elements instead of shells in modeling buckling (apart from the lack of remote loading availability which the shell edges in CalculiX, of course…)

Ok, now I am more curious to see the results of shell elements (with a single layer), and eager to share our findings…

[Calc_em]

@fethio

So, S4R expands to a single layer of C3D8R, which suggests that it will be advisable to switch to a volume model and create more than one layer in the thickness direction, right?

Yes, but you could also use a multilayer (same material) composite shell section supported for S6 and S8R elements, or simply try with those second-order elements right away.

[fethio]

Going back to basics, we built a cylinder without holes and trying various element types and formulations, mesh densities, boundary conditions in modeling buckling of a thin-walled long cylindrical vessel under axial compression.

The reference we are using to calculate critical buckling load is NASA/SP-8007-2020/REV 2 (Though we are not showing any calculations in this post.)

We figured that the previous model in this thread had a fixed-free boundary condition. We modified it to reflect the simply supported case to conform with the above reference. (We assumed that the simply-support condition corresponded to restraining the radial degrees-of-freedom of the nodes sets lying along the TOP as well as the BOTTOM egdes.)

We seek for the buckling modes. A small axial compressive load is used as the perturbation load. TRANSFORM card is used to align DOFs of the top and bottom edge nodes with the cylindrical r,theta,z coordinates. BOUNDARY card is used to restrain the DOF-1 (r) of both node sets, while DOF-3 (z) is additionally restrained at the bottom.

*Heading
linear buckling analysis, skirt without openings and holes under load, shell mesh lc 128, quad 4.

*INCLUDE, INPUT=noholes_S4_128.inp
*TRANSFORM, NSET=skir_sless_top_LOAD, TYPE=C
0.0, 0.0, 0.0, 0.0, 0.0, 1.0
*TRANSFORM, NSET=base_plate_center_YFIX, TYPE=C
0.0, 0.0, 0.0, 0.0, 0.0, 1.0
**TRANSFORM, NSET=skir_carbo_und, TYPE=C
**0.0, 0.0, 0.0, 0.0, 0.0, 1.0

*BOUNDARY
base_plate_center_YFIX,1,1
base_plate_center_YFIX,3,3
skir_sless_top_LOAD,1,1

*MATERIAL,NAME=STEEL
*ELASTIC
2.09000000e+05, 2.80000000e-01
*DENSITY
7.8E-9

*ORIENTATION, NAME=CYL_CSYS, SYSTEM=C
0.0, 0.0, 0.0, 0.0, 0.0, 1.0

*SHELL SECTION, ELSET=skir_sless, MATERIAL=STEEL
25
*SHELL SECTION, ELSET=skir_carbo_mid, MATERIAL=STEEL
25
*SHELL SECTION, ELSET=skir_carbo_und, MATERIAL=STEEL
25

** step 1, linear buckling using preloaded state
*STEP, PERTURBATION
*BUCKLE
10
*CLOAD
skir_sless_top_LOAD,3,2337.5
*NODE FILE, OUTPUT=2D
U
*END STEP

Based on a mesh of S4 to start with, out of the 10 buckling modes solved for, all modes exhibit a homogeneously radially pulsating cylinder. For one thing, the buckling mode does not seem to reflect the imposed boundary conditions. I would expect the radial displacements to be zero at the top and bottom edges but they are not. Secondly, I expected to see a variety of modes while all modes converged to the identical mode.

https://limewire.com/d/UcAUN#W8VQMgkSWx

Edit: My gut feeling tells me that the transformation of coordinates may not be working as expected for the shell type, due to the extra rotational DOFs per node. I would appreciate any insights on this.

[Calc_em]

With the current BCs, the cylindrical shell is free to move in the tangential direction (rotate around its axis). Make sure that the model isn’t underconstrained in the buckling analysis (frequency analysis can show it well).

[fethio]

Thank you for pointing out the underconstrained condition @Calc_em. As @Disla suggested we went back to working with cylinders without holes under axial compressive action (with analytical formulae to calculate critical buckling modes).
From NASA/SP-8007-2020/REV 2 :
4.1.1.1 Axial Compression (p.22)

	      % input model parameters (all dimensions in m's)
Len = 10.5     % cylinder height
Rad = 1.6275   % cylinder shell mid-shell radius
Dia = 2*Rad    % diameter
thk = .025     % cylinder shell thickness
fa1 = Rad/thk  % frequently used ratio

               % input material properties
E_m = 2.09E11  % Elastic modulus, in Pa
nu_ = 0.28    % Poisson's ratio		      
fa2 = (1-nu_^2) % frequently used factor

               % buckling mode geometry factors
m_ = 10         % the number of axial half waves in the mode
n_ = 10         % the number of circumferential full waves
gam = 1       % buckling knock-down factor

               % structural properties
              % is the wall flexural stiffness per unit width (Eq.3)
D_ = E_m*thk^3/(12*fa2)
              % curvature parameter (Eq.4)
Z_ = Len^2/(Rad*thk)*sqrt(fa2)
              % buckling aspect ratio (Eq.5)
bet = n_*Len/(m_*pi*Rad)
fa3 = m_^2 * (1+bet^2)^2 % frequently used factor
              % buckling coefficient  (Eq.2)
k_x = fa3 + (12/pi^4)*(gam*Z_)^2/fa3

              %  buckling line load of a simply supported cylinder 
              % under axial compression
              % based on Donnell’s shell theory (Eq.1)
N_x = k_x*(pi/Len)^2*Dia
             % total buckling load
F_x = (2*pi*Rad)*N_x

              % cross check w/ reduced form
             % critical buckling stress (Eq.7)
s_x = gam*E_m/sqrt(3*fa2)*1/fa1
             % transverse cross-sectional area
ar_ = pi*Dia*thk
             % calculation of percentage difference
             % between the two 
abs(s_x*ar_ - F_x)/F_x

I cannot figure out the reason for the HUGE difference between critical buckling loads found from 1) the buckling line load Eq.5, and 2) critical buckling stress formula, Eq.7. It turns out the there’s a difference of nearly four orders of magnitude (one is roughly 5000 times the other)

FEA solution that I post below favors the second (reduced, Eq.7) formula. The model subjected to axial displacement of 100 mm’s using S3, S4R, S6, S8, and S8R element types. The boundary conditions were fixed at the base and free at the top. (We are working on the simply supported). The convergence characteristics of the elements (number of iterations, final displacements, and corresponding reaction forces) are

| Element Type | Increment | Displacement |           RF3 |
|--------------+-----------+--------------+---------------|
| S3           |        21 |         97.4 | -4.805185E+08 |
| S4           |        14 |         49.7 | -2.452951E+08 |
| S4R          |        21 |         39.0 | -1.921596E+08 |
| S6           |        18 |         51.1 | -2.543565E+08 |
| S8           |        11 |         50.7 | -2.523001E+08 |
| S8R          |        20 |         46.9 | -2.332738E+08 |

The buckling load calculated from Eq.7 is 4.94E+08 N. I am missing something in the calculation based on Donnell’s shell theory, or maybe a typo in the NASA manual !!

By the way, all above solution types diverge at the verge of a possible buckling mode, but the geometry remains isomorphic with respect to the reference (perfect cylinder) geometry, i.e. there are no signs of any buckling deformation until the end of all analyses.

N.B. The Donnell’s formula produces a buckling load that is similar to Eq.7 if we use m=n=0.01.

[Disla]

Aren’t you mixing units?

Eq.7 delivers critical buckling stress (MPa) not (N)

[fethio]

I am multiplying the outcome of Eq.7 (critical buckling stress) with the tube’s cross-sectional area to convert it to the critical buckling load, but still a difference of about four orders exist…

[Calc_em]

@fethio

N_x = k_x*(pi/Len)^2*Dia

Isn’t D in this formula supposed to be the stiffness and not the diameter ?

[fethio]

Of course it is, thanks @Calc_em . I should have looked more carefully.

now with m=1/2 (a quarter of a full axial wave) and n=2 (two full circumferential waves) the two formulas are in agreement as close as 3%.

And the first mode shape using S4 element looks like this:

[Disla]

Hi. Donell formula is straight forward. It’s validity range assumes thin shell which means you can expect increasing discrepancy with your fea with increasing thickness.

To nail the right set up it will be more comvinient to use 5mm which should provide better agreement.

m and n values are strongly degenerated and valid for comparison to find out if your solution makes sense.They are integer values. Displacement field should be symmetric and yours is not meaning you BC are not right yet or you are usng different BC than the ones Donell formula..

This are my results for 5mm and 10mm. S8R and Kinematic Couplings . Simply supported allows rotation of the shell.

[fethio]

Hi @Disla , just to understand your model setup, are you applying a preload such that it produces an axial stress of 1 MPa?

And concerning the boundary conditions, I understand that having the z-axis aligned with the longitudinal axis, you fixing all the transverse (x,y) translational degrees-of-freedom on the upper and lower edges, while the remaining three rotations are free. In our model we restricted to restrain only the radial displacements (after transforming the edge node set degrees-of-freedom to cylindrical system).

[Calc_em]

@fethio

In our model we restricted to restrain only the radial displacements (after transforming the edge node set degrees-of-freedom to cylindrical system).

As discussed before, it’s not suffcient and you need some additional constraints (even limited to individual nodes) to prevent the remaining rigid body motions. If only symmetric mode shapes are needed, you could utilize symmetry, but it doesn’t work well with shells in CalculiX. Did you try with solid element layers ?

[fethio]

apart from the transverse translational restraints we are applying an axial translation restraint on one of the edge node set and an additional circumferential restraint on a node that lies on either one of the edge node sets.

By the way, looking at the previous discussion threads about applying remote loading to shell edges, does kinematic coupling work on shell edges? I recall there were some problems about making this work.

[Disla]

@fethio

Hi @Disla , just to understand your model setup, are you applying a preload such that it produces an axial stress of 1 MPa?

Right. You can also consider force applied to the kinematic ref NODE. Your Buckling load will be delivered in terms of Fcr (N)

@fethio

on a node that lies on either one of the edge node sets

That breaks the simmetry of the displacement field. Try with two in oposite sides to get a cleaner solution.

@fethio

does kinematic coupling work on shell edges

In my experience yes as far as it is defined in terms of a SURFACE TYPE=NODE.

*SURFACE,NAME=NSet_KINEMATIC,TYPE=NODE
NSet_KINEMATIC
*COUPLING,REF NODE=XXX,SURFACE=NSet_KINEMATIC,CONSTRAINT NAME=RootCouple
*KINEMATIC
1,3

Don’t forget to run a frequency analisys before Buckling. Kinematic needs a special care to avoid underconstraining the model. Keep this image in your brain when using KInematic.

Prepomax has implement it in it’s last version. Would be a good oportunity to test it althought a saw something strange in the example. I will post there to comment.

[Calc_em]

@fethio

By the way, looking at the previous discussion threads about applying remote loading to shell edges, does kinematic coupling work on shell edges?

Couplings and rigid bodies may not work well with shells in nonlinear analyses: Rigid-body constraint convergance problems - #25 by dhondt

In linear ones, they are usually safe.

In Abaqus, it’s usually better to use distributing couplings with shells.

Also see this thread: Different coupling constraints and their limitations

Element-based distributing coupling is more difficult to define, but may help avoid the aforementioned convergence (and potentially other) issues.

[fethio]

Thank you for the suggestions. I am not using the *COUPLING or *RIGID BODY cards here.

We obtained a solution that is in agreement with our reference formula (derived from Donnell’s shell formulation). The model uses S4 elements, simply supported end conditions (FIXED x,y,z translations at the BOTTOM edge and FIXED x,y translations at the TOP edge). The TOP edge is loaded so as to give an axial compressive stress of 1 MPa only. Thus the buckling factor is directly interpreted as the critical buckling stress.

In our solution the 4th linear buckling mode gives the picture below.

The critical buckling stress reads 1929 MPa. From the picture the number full circumferential waves is n=5 and the number of half axial waves is m=4. Donnell’s formula using m=4 and n=5 yields a buckling load of 4.9362e+08 N. Eq.7 gives 4.9360e+08, that is a perfect agreement.

@Disla

@Disla , Comparing with your solution, I suspect that increased thickness in my model (25 mm vs 5 or 10 mm in yours) causes smaller m and n values for the fundamental buckling mode?

[Disla]

Yes, you are absolutely right. The smaller the thickness the larger the m,n numbers.

Looks like you need to increase the accurcy a little (or request more number of modes) to get a better mode shape.

Accuracy goes as 1/n^2 meaning the smaller the thickness (Larger number of waves) the better agreement with Donell. Do not forget the larger the number of waves , the samer should be your mesh to capture properly the pattern.

¿Why are you showing/interested in the 4th mode?