Shell Elements for Stress Analysis

1. Introduction

The shell elements are designed for static and dynamic, linear and nonlinear analysis of thin- and thick-walled plate and shell structures with isotropic, orthotropic and laminated, linear and nonlinear materials. They are based on Reissner-Mindlin shell theory and have 5 or 6 degrees of freedom per node. To account for geometric nonlinearity, a Total-Lagrangian approach is employed. The MITC (Mixed Interpolation of Tensorial Components) formulation[BatheMITC] improves the out-of-plane behaviour by interpolating the through-the-thickness term of the strain tensor differently. This is a better way to prevent out-of-plane locking than selective under-integration. The rotation of the director is computed by making use of finite rotations to allow for computing the director of the deformed configuration directly by an algebraic expression of the rotational degree of freedoms, i.e without numeric integration.

For laminates involving linear materials for all layers, the shell elements automatically perform the through-the-thickness integration a priori in analytical fashion. Hence, the number of layers is irrelevant to the computation time of the first and second variations. The benefit of this optimization is a considerable speed-up for analysis of composite structures. Note that this optimization does not alter the element response.

2. Element Types

The following element types are available:

Table 17. Shell Elements for Stress Analysis

Element type Remarks
T3.S.MITC  
T6.S.MITC  
Q4.S.MITC  
Q4.S.MITC.E4 Like the Q4.S.MITC element, but with 4 incompatible modes for in-plane bending.
Q8.S.MITC  
Q9.S.MITC  


The following points are worth mentioning:

  • The Q9.S.MITC element is more effective (i.e. it has a better ratio of accuracy vs. computational effort) than the other B2000++ shell elements, even for highly-distorted meshes. This element should be considered whenever a new FE model is being made. The Q8.S.MITC element is intended for FE models that were imported from other FE codes that do not support nine-node quadrilateral elements.

  • Due to the incompatible modes, the Q4.S.MITC.E4 element does not pass the patch tests (i.e. cannot exactly represent constant strain fields), and it may exhibit instabilities in buckling analysis and in nonlinear analysis. Elements of this type should not be significantly distorted. While the in-plane bending response of this element is better than the Q4.S.MITC element, the latter gives more predictable results, especially concerning the gradients.

  • For the the quadratic MITC elements (T6.S.MITC, Q8.S.MITC, Q9.S.MITC), mid-edge and mid-face nodes should be located at the middle of the edges and centre of the face, respectively. Otherwise, the Jacobian will contain higher-order terms, and the elements will not work with full accuracy. However, for moderately curved geometries such as cylindrical panels, the error is small enough that the accuracy is still very good.

  • The T6.S.MITC element fails for certain double-curved geometries such as hyperboloid shells. Further, for complex loading conditions, this element may exhibit instabilities which affect the accuracy of the linearised buckling load (i.e. the predicted buckling load will be too low for coarse meshes).

3. Required element attributes

mid m

Specifies the element material number m. The elements can process materials of the following types:

nthickness t1 t2 t3...

Specifies the element's thickness at each of the element's nodes. This option allows to specify varying element thickness. The element thickness is interpolated using the element's in-plane shape functions. Not required for laminate materials.

thickness t

Specifies a constant element thickness t. Not required for laminate materials.

4. Optional element attributes

eccent e

Defines the element eccentricity e, which is a constant offset of the shell volume w.r.t. the shell reference surface. The same units as for the coordinates and thicknesses are used. Default is 0. With e = 1 2 t , where t is the shell thickness, the shell bottom surface coincides with the shell reference surface. With e = 1 2 t , the shell top surface coincides with the shell reference surface.

group eid

Defines the element group number gid (a non-negative integer number). The default group number is 0. The same definition will be used for all elements defined hereafter, until a new group option is encountered or until the eltype command is specified.

morientation default

Specifies the default material orientation which is obtained by projecting the branch-global x-direction onto the shell surface and constructing an orthogonal reference frame (this will not work for elements whose shell surface normal is aligned with the x-direction).

morientation ... end

Specifies the material orientation. This can be done in one of several ways.

The first way is to specify a base with two vectors, from which an orthogonal reference frame is constructed (the base specified base vectors do not need to be orthogonal, but they must not be colinear). This reference frame can be rotated about one of its axes. One of the axes is then projected onto the element reference surface. Using the shell surface normal, an orthogonal reference frame is constructed. Finally, this reference frame can be rotated about the element surface normal:

morientation
  base u1 u2 u3  v1 v2 v3
  [rotate axis X|Y|Z angle a]
  project axis X|Y|Z
  [rotate angle a]
end

The second way specifies a transformation, from which an orthogonal reference frame for the coordinates of the element integration point is calculated. This reference frame can be rotated about one of its axes. One of the axes is then projected onto the element reference surface. Using the shell surface normal, an orthogonal reference frame is constructed. Finally this reference frame can be rotated about the element surface normal:

morientation
  transformation id
  [rotate axis X|Y|Z angle a]
  project axis X|Y|Z
  [rotate angle a]
end

The third way specifies a vector which is projected onto the shell reference surface. Using the element surface normal, an orthogonal reference frame is constructed, whose Z-axis is aligned with the element surface normal. This reference frame can then be rotated about the element surface normal:

morientation
  vector u1 u2 u3
  [rotate angle a]
end

Yet another possibility is to use the element-local reference frame as calculated at the element integration point (the Z-axis is the shell surface normal). This reference frame can be rotated about the shell surface normal:

morientation
  element
  [rotate angle a]
end

The rotation angles are specified in degrees. When using laminate materials, a final rotation about the ply angle performed automatically. For the projection, it is necessary that the projected axis or vector is not colinear with the shell surface normal.

To visualize with baspl++ the material orientations as calculated at the element integration points, an analysis (for example linear) needs to be run, and with the option gradients enabled. This will write the dataset MBASE_IP to the database which can be visualized by means of the following script:

if len(sys.argv) != 2:
    print 'usage: baspl++ %s database' % sys.argv[0]
    sys.exit(1)

m = Model(sys.argv[1])
p = NPart(m)
p.edge.show = True
p.elements.extract = 1
f = Field(m, 'MBASE_IP', case_index=0, cycle_index=0)
p.sampling.field = f
neccent e1 e2 e3 ...

Defines the nodal eccentricities e1, which are the offset of the shell volume w.r.t. the shell reference surface, at the nodal position(s). This command allows to specify varying eccentricity. The element eccentricity is interpolated using the element's in-plane shape functions. Otherwise, the definition is similar to eccent.

non_structural_mass v

Defines the non-structural mass per unit area. Default is 0. The same definition will be used for all elements defined hereafter, until a new non_structural_mass option is encountered or until the eltype command is specified.

5. Optional case attributes

drills f

Defines the drill stiffness factor for shell element nodes with 6 degrees-of-freedom and whose in-plane rotation (drill) will be stabilized. Default is 1e-8. The effective drill stiffness is calculated using the averaged in-plane shear modulus and the element configuration and thickness, and is multiplied with f. Hence, a value of 0 for f means no drill stabilization.

Compared to the autospc option in conjunction with the MUMPS sparse linear solver (see Sparse Linear Solvers), drills applies only to the 6th degrees-of-freedom of shell elements that may need stabilization.

6. 5 or 6 degrees-of-freedom

A well-known problem with shell elements is how intersections or junctions (such as the connecting edge between a panel and a stiffeners) are treated. Pure 5 degrees-of-freedom elements are not able to take into account the drill moment and the rotation in z-direction. This violation of the compatibility requirement results in numerical problems that often prevent a successful simulation. Elements with an additional degree-of-freedom are able to maintain compatibility at junctions. However, in places where there is no such junction, 6 degrees-of-freedom Elements must introduce an artificial drill stiffness term to prevent the second variation from becoming singular.

The MITC elements largely eliminate the need for artificial drill stiffness by allowing for each node to have either 5 or 6 degrees-of-freedom. The B2000++ input processor b2ip++ determines automatically which shell element nodes have 5 degrees-of-freedom, which ones have 6 degrees-of-freedom (see the MDL commands detect_shell_intersections and shell_intersection_angle), and which ones of the latter will be drill-stabilized. This is done such that junction nodes and nodes where boundary conditions or linear multipoint constraints on one or several rotational degrees-of-freedom are defined have 6 degrees-of-freedom, while the other nodes have 5 degrees-of-freedom. Nodes that connect to beam elements also have 6 degrees-of-freedom. The shell element nodes that are marked to have 6 degrees-of-freedom can be visualized in the post-processor baspl++ by selecting the node type "19":

if len(sys.argv) != 2:
    print 'usage: baspl++ %s database' % sys.argv[0]
    sys.exit(1)

m = Model(sys.argv[1])
p1 = NPart(m)
p1.nodes.ntypes = 19
p1.nodes.extract = True
p2 = NPart(m)
p2.edge.show = True
p2.opacity = 0.1
p2.elements.extract = True

Nodes whose drill rotation will be stabilized are restricted to the smallest possible set. For instance, nodes at intersections do not need stabilization, because the adjacent shell element(s) will provide the required stiffness. Nodes that are connected to a beam element do not need stabilization, since the beam element will provide the required stiffness. For the nodes that will be stabilized, the input processor b2ip++ will create the nodeset "b2ip_drills" which can be visualized with the following script:

if len(sys.argv) != 2:
    print 'usage: baspl++ %s database' % sys.argv[0]
    sys.exit(1)

m = Model(sys.argv[1])
p1 = NPart(m)
p1.nodes.nodesets = 'b2ip_drills'
p1.nodes.extract = True
p2 = NPart(m)
p2.edge.show = True
p2.opacity = 0.1
p2.elements.extract = True

This approach ensures that most of the numerical problems traditionally associated with the 6th degree-of-freedom are resolved. In addition, for the nodes with 5 degrees-of-freedom, node normal vectors are computed by the input processor. Having continuity of the shell normal between elements for proper handling of shell intersections and correctly transferring twisting moments in curved shell models [Hoff95], ensures that adjacent elements are truly compatible along the edges and is a prerequisite for solving problems like the Raasch Challenge [Knight96].

This works even when a part of the shell elements is inverted such that the shell element normal points in the opposite direction. However, opposite shell element normals may be due to modelling errors and may cause problems in conjunction with non-symmetric laminates or traction loads. For this reason, the input processor b2ip++ emits a warning when opposite shell element normals are detected, and will create the nodeset "b2ip_opposite" which can be visualized with the following script:

if len(sys.argv) != 2:
    print 'usage: baspl++ %s database' % sys.argv[0]
    sys.exit(1)

m = Model(sys.argv[1])
p1 = NPart(m)
p1.nodes.nodesets = 'b2ip_opposite'
p1.nodes.extract = True
p2 = NPart(m)
p2.edge.show = True
p2.opacity = 0.1
p2.elements.extract = True

7. Stresses and Strains

Stresses, strains, and any failure criteria which are stored on database, are computed at the integration points ("Gauss" points) which are used by the element to compute the first and second variations. The stresses and strains which are stored on database are all expressed in the branch-local coordinate system. For linear analysis, the strains and stresses stored on the database are the "engineering" stresses and strains. For nonlinear analysis, the strain is the Green-Lagrange strain, while the stress is the Cauchy stress.

The in-plane integration rule depends on the element type. For the through-the-thickness integration, a 3-point Gauss-Legendre-Lobatto rule is used (for laminates, this applies to each layer). This rule ensures that gradients are evaluated at the surfaces and at the layer boundaries, thus, where the maximum strains and stresses occur.