Shell Elements


What is a Shell Element?

Shell elements are 4- to 8-node isoparametric quadrilaterals or 3- to 6-node triangular elements in any 3-D orientation. The 4-node elements require a much finer mesh than the 8-node elements to give convergent displacements and stresses in models involving out-of-plane bending. Figure 1 shows some typical shell elements.

 

The General and Co-rotational shell element is formulated based on works by Ahmad, Iron and Zienkiewicz and later refined by Bathe and Balourchi. It can be applied to model both thick and thin shell problems. Also, the geometry of a doubly curved shell with variable thickness can be accurately described using this shell element.

 

Figure 1: Sample 3-D Shell Element
 

Figure 2: Typical Shell Elements

 

The Thin shell element is based on thin plate theory.  The bending behavior of the element is based on a discrete Kirchoff approach to plate bending using Batoz's interpolation functions.  This formulation satisfies the Kirchoff constraints along the boundary and provides linear variation of curvature through the element.  The membrane behavior of the element is based on the Allman triangle which is derived from the Linear Strain Triangular (LST) element.  A general curved surface is approximated by this element as a set of facets formed by the planes defined by the three nodes of each element.  For these reason a well-refined mesh is necessary.

 

The element geometry is described by the nodal point coordinates. Each shell element node has 5 degrees of freedom (DOF) - three translations and two rotations. The translational DOF are in the global Cartesian coordinate system. The rotations are about two orthogonal axes on the shell surface defined at each node. The rotational boundary condition restraints and applied moments also refer to this nodal rotational system. The two rotational axes (V1 and V2) are usually automatically determined by the processor and you do not have to specifically orient them.

 

The rotational R1 and R2 directions at each node are determined by the cross product of the normal vector Vn and a guiding vector Vg as follows. See Figure 3 (V1 and V2 are the unit vectors along the R1 and R2 rotational direction, respectively):

R1 = V2 x Vn

R2 = Vn x Vg

 

The default guiding vector Vg is (1,0,0). In other words, the R1 axis is actually the projection of the guiding vector Vg onto the shell surface, and the R2 axis lies on this shell surface but orthogonal to R1 and the normal vector Vn. Figure 3 shows the orientation of R1 and R2 using the default guiding vector (1,0,0).

 

You can enter two guiding vectors for each shell element part . If the guiding vector has zero projection onto the shell plane, the second guiding vector is used. You can also specify a guiding vector for a specific node.

 

The processor determines the normal direction of each shell element node by taking the average of different normal vectors from different elements connecting to the node point. The element normals are computed from the element corner nodes data. From this normal vector, two rotational axes are then determined. You can also specify the normal vector direction directly and suppress the processor's calculations. The two rotational vectors can also be oriented in a user-specified direction such that skewed rotational restraints or moments in 3-D space can be easily defined.

 

When multiple shell panels intersect at an area (usually at a line of intersection), it is difficult to specify the rotational DOF at the area. In this case, the nodes at the intersecting area should be designated as 6-DOF nodes.
 


Figure 3: The Rotational Degrees of Freedom Determined Using Vg

 

Six DOF nodes will have globally oriented rotational DOF instead of the unique R1 and R2 DOF at each node.

 

Three-dimensional shell elements are Type 26 elements and are 4- to 8-node isoparametric quadrilaterals or 3- to 6-node triangular elements in any 3-D orientation. They can be applied to models for both thick and thin shell problems. Also, the geometry of a doubly curved shell with variable thickness can be accurately described using this shell element. Theoretical aspects of this element formulation can be found in Bathe and Balourchi.

 

The element geometry is described by the nodal point coordinates. Each shell element node has five degrees of freedom (DOF) - three translations and two rotations. The shell elements' rotational DOF are defined by the shell guiding vectors. The projection of the guiding vector, Vg, onto the shell element will be the first rotational axis (R1). The second rotational axis is orthogonal to both the R1 and the shell normal vector (Vn), using the rule R2 = Vn x Vg. The default guiding vector is (1, 0, 0). If the guiding vector is parallel to the shell normal vector, Vn, a second guiding vector [(0, 1, 0) by default] will be used to define the R1 direction.


Shell Element Parameters

There are three types of shell elements available that provide different options. The difference between these types is explained in the What is a Shell Element? section. The shell element type can be selected in the "Element Formulation" drop-down box.

 

First, you must specify the material model for this part in the "Material Model" drop-down box. The available shell material models are:

Tip

Using a higher integration order (set on the "Advanced" tab) for hyperelastic materials can help with convergence when the part experiences large deformation.

 

Specify the thickness of the shell elements in the "Thickness" field.The element is considered to be drawn at the midplane of the shell element. You must enter a value for the thickness to run the analysis. If the "Composite" option is selected in the "Material Model" drop-down box, this field will not be available. In this case, the sum of the thicknesses of the laminae defined in the "Composite" tab will be displayed in this field.

 

If you want the shell elements in this part to have the midside nodes activated, select the "Included" option in the "Midside Nodes" drop-down box. If this option is selected, the shell elements will have additional nodes defined at the midpoints of each edge. (For meshes of CAD solid models, the midside nodes follow the original curvature of the CAD surface, depending on the option selected before creating the mesh. For hand-built models and CAD model meshes that are altered, the midside node is located at the midpoint between the corner nodes.) This will change a 4-node shell element into a 8-node shell element. An element with midside nodes will result in more accurately calculated gradients. Elements with midside nodes increase processing time. If the mesh is sufficiently small, then midside nodes may not provide any significant increase in accuracy.

 

Tip

The displacement at midside nodes are always output. The stress and strain at midside nodes are output only if the user activates the option to output these results before running the analysis. The option is located under the "Analysis: Parameters...: Advanced" dialog on the "Output" tab. (See the page "Setting Up and Performing the Analysis: Nonlinear: Analysis Parameters: Advanced Settings: Controlling the Output Files" for details.)

 

Use the "Analysis Type" drop-down to set the type of displacement that is expected. "Small Displacement" is appropriate for parts that experience no motion and only small strains and will ignore nonlinear geometric effects that result from large deformation. "Large Displacement" is appropriate for parts that experience motion and/or large strains. (Depending on the type of shell element chosen, there may be other options on the "Analysis Formulation" on the "Advanced" tab that are affected by the "Analysis Type" drop-down.)

 

Tip

Use the "Tools: Options: Analysis" tab and set the "Use large displacement as default for nonlinear analyses" option to control whether the "Analysis Type" defaults to small or large displacement.


 

Using the Composite Material Model

If the shell elements are using a composite material model, you can select a failure criteria using the "Failure Criterion" drop-down box in the "Composite" tab of the "Element Definition" dialog. The available options and the corresponding equations used are:

(1)

 

where

 

Xt = Axial or longitudinal strength in tension (>0)

Xc = Axial or longitudinal strength in compression (>0)

Yt = Transverse strength in tension (>0)

Yc = Transverse strength in compression (>0)

S = Shear strength

F12 = the stress interaction value entered in the material properties. If the value is 0, it will default to

s1 = Stress in principal material 1- direction

s2 = Stress in principal material 2- direction

t12 = Shear stress in principal material 1-2 plane

 

If equation (1) is not satisfied, the material has failed.

 

(2)

where R is the strength/stress ratio.

Since each combination of stress components in (1) reaches its maximum when the left-hand side reaches unity, (2) can be substituted into (1) and yield:

Solving this equation for R yields:

where:

with the signs of the stresses sij reversed.

(3)

 

where

s1 = Calculated stress in 1 direction

s2 = Calculated stress in 2 direction

t12 = Calculated shear stress

Xc = Allowable compressive stress in 1 direction (>0)

Yc = Allowable compressive stress in 2 direction (>0)

Xt = Allowable tensile stress in 1 direction (>0)

Yt = Allowable tensile stress in 2 direction (>0)

S = Allowable shear stress (>0)

 

If the three conditions of equation (3) are not satisfied, the material has failed.

(4)

 

where

e1 = Calculated strain in 1 direction

e1 = Calculated strain in 2 direction

g12 = Calculated shear strain

T1c = Allowable compressive strain in 1 direction (>0)

T2c = Allowable compressive strain in 2 direction (>0)

T1t = Allowable tensile strain in 1 direction (>0)

T2t = Allowable tensile strain in 2 direction (>0)

S = Allowable shear strain (>0)

 

If the conditions of equation (4) are not satisfied, the material has failed.

 

The "Composite Laminate Stacking Sequence" table can be used to define the layers of the composite. Layers are oriented and ordered as described below in "Controlling the Orientation of Shell Elements". The "Thickness" column must be defined for each layer. The "Orientation Angle" column gives the angle a between the material axes (defined on the "General" tab) and the fiber or layer axes. The material properties are entered according to the layer axes. To define the material properties for a layer, click in the "Material" column. A dialog will appear that will allow to select an existing material or "Add" a new material. (Also see the page "Setting Up and Performing the Analysis: Nonlinear: Material Properties: Composite Material Properties" for details.)

 

Tip

Use the "Copy Row" button to copy existing input to create repeating series. Either enter a single lamina to copy, or a range of lamina separated by a dash "-", such as 4-8. For the "Insert Before Row" field, specify any row number after the last entry in order to insert the copy at the end of the stacking sequence.


 

Defining Thermal Properties of Shell Elements

If a shell element part is using a thermoelastic or  temperature-dependent material model, you must specify a value in the "Stress free reference temperature" field in the "Thermal" tab of the "Element Definition" dialog.  This value is used as the reference temperature to calculate element-based loads associated with constraint of thermal growth using bilinear interpolation of the nodal temperatures.

 

If a shell element part is using a viscoelastic material model and you want creep effects to be analyzed, select the "Power-law" option in the "Creep law" drop-down box.  This will cause the power law creep function to be used to calculate the creep effects during the analysis.  If you want the creep calculations to be calculated on evenly sized divisions of the time step, select the "Fixed substeps" option in the "Time integration method" drop-down box.   If you want the creep calculations to be calculated on variable sized divisions of the time step, select the "Flexible substeps" option.  Specify the temperature at which no thermal stress exists in the "Stress free reference temperature" field.  If you are performing an analysis with non-cyclical loading, select the "Effective" option in the "Creep strain definition" drop-down box.  If you are performing an analysis with cyclical loading, select the "Modified" option.  

 

During the analysis, the creep calculations will be performed as iterations in substeps of each time step.  You can control how many substeps are allowed in a single time step in the "Maximum number of substeps" field.  You can also specify how many iterations can be performed in a single substep in the "Maximum number of iterations in a substep" field.  After each substep iteration, the creep stress and strain will be compared to the previous iteration.  If the value is not within the tolerances specified in the "Creep strain calculation tolerance" and "Creep stress calculation tolerance" fields, another iteration will be required.  If you want to use a fully explicit method for the itme-0integration scheme, type 0.0 in the "Time integration parameter" field.  If you want to use a fully implicit method, type 1.0 in the "Time integration parameter" field.


Controlling the Orientation of Shell Elements

There are two types of axes that the user can control. One type is the axis normal to the shell elements. The main reason to specify the normal direction is for applying pressure to the element or specifying surface-to-surface contact. For composites, the normal direction controls which side of the element layer 1 is on. The second type of axes is the element's in-plane axes. This is useful when working with an orthotropic material model or you desire to get results in the element coordinate system.

 

In addition to these two types, there are three axis systems which the user may need to set:

  1. The element axes.

  2. The material axes, which is used for orthotropic and composite material models.

  3. The laminate axes, which is used with composite material models.

Element Axes:

An element normal point is used to control the orientation of the element's normal axis (+3), or which side of the element is the top side (+3 side) and the bottom side (-3 side). The normal direction is determined by specifying a point in space using the "X Coordinate", "Y Coordinate" and "Z Coordinate" fields in the "Element Normal" section of the "Orientation" tab. See Figure 4. A positive normal pressure will be applied normal to the shell elements the direction of the +3 axis and thus will point away from the element normal point.

 

(X,Y,Z)

Element Normal

Coordinate

 

 

bottom side

+3 axis

top side

Figure 4: Determining the Element Normal
The edge-on view of the shell element is shown.

 

For a general FEA analysis, you can ignore the element's in-plane orientation, axes 1 and 2. The ability to orient element's in-plane axes is useful for elements with orthotropic material models and composite elements. The orientation of the 1 and 2 axes is done in the "Orientation" tab of the "Element Definition" dialog. The "Method" drop-down box contains three options that can be used to specify which side of the element will be the ij side. The element axis 1 will be parallel to the ij side of the element. (Element axis 2 forms a right-hand system with axes 1 and 3.)

Material Axis:

The material properties are entered in terms of the material axes for orthotropic materials, and the material axes can be defined from the element axes (or from a global reference) using the "Material Axis Direction" section of the "General" tab on the Element Definition dialog. (Material properties for composite materials are entered in terms of the layer axes, as described below.)

 

In all cases, the material axis 3 will be perpendicular to element and in the same direction as the element axis 3.

 

If the Co-rotational shell elements are using a orthotropic or composite material model, the in-plane material axes 1 and 2 can be defined as follows:

Regardless of the method used to orient the Co-rotational shell element material axes, you can use the "In-Plane Rotation Angle" field to rotate the material axes by a specified angle about axis 3. The rotations follow the right-hand rule about axis 3.

 

If the General shell elements are using an orthotropic material model, the in-plane material axes 1 and 2 can be oriented using two methods. The first method is to specify a vector parallel to material axis 1 by using the "X Direction", "Y Direction" and "Z Direction" fields. Material axis 2 is in the plane of the element and forms a right-hand system with axes 1 and 3. The second method is to define an angle in the "Material Axis Rotation Angle" field that will rotate material axis 1 about the element 3 axis using the right-hand rule (see Figure 5). This method will only be used if the values in the "X Direction", "Y Direction" and "Z Direction" fields are 0.

 


Figure 5: Rotation Angle Definition of Material Axis Direction

Layer Axes:

For composites, each layer of the laminate has a set of local axes. The layer 1-axis is the axis along the fiber of each individual layer. The layer 2-axis is the axis perpendicular to the fiber of each individual layer and in the plane of the element. The layer 3-axis is perpendicular to the element and therefore parallel to the local element axis c. (To avoid confusion between the layer axes and the element axes, the element axes are often referred to as a-b-c when working with composites, and axes 1-2-3 define the orientation of the fibers.) See Figures 6 and 7.

 

Figure 6: Laminae Stacking Order

The user enters the angle a for each layer using the "Orientation Angle" of the "Composite Laminate Stacking Sequence" table. Axis 3 is controlled by the element normal coordinate.

 

Figure 7: Alternate Depiction of Element, Material, and Layer Axes

The element axes a-b-c are based off of the element normal coordinate and the ij side of the element. The material axes x-y-z are rotated an angle b from the element axes. The element and material axes are the same for all layers in the stack. The layer or fiber axes 1-2-3 are rotated an angle a from the material axes and can be different angles for each lamina in the stack.


Advanced Shell Element Parameters

Select the formulation method that you want to use for the shell elements in the "Analysis Formulation" drop-down box in the "Advanced" tab.

If the "Allow for overlapping elements" checkbox is activated, overlapping elements will be allowed to be created when the lines are decoded into elements. Overlapping may be necessary when modeling elements. This is especially true for problems confined to planar motion.

 

If you want the stress results for each element to be written to the text log file at each time step during the analysis, activate the "Detailed stress output" checkbox.  This may result in large amounts of output.

 

If one of the von Mises material models has been selected, you can choose to have the current material state (elastic or plastic), current yield stress limit, current equivalent stress limit and equivalent plastic strain output at corner nodes and/or integration points at every time step.  This is done by selecting the appropriate option in the "Additional output" drop-down box.

 

If you are using General shell elements, select the integration order that will be used in this part in the "Integration Order" drop-down box.  For rectangular shaped elements, select the "2nd Order" option.  For moderately distorted elements, select the "3rd Order" option.  For extremely distorted elements, select the "4th Order" option.  The computation time for element stiffness formulation increases as the third power of the integration order.  Consequently, the lowest integration order which produces acceptable results should be used to reduce processing time.  

 

The reduced integration scheme can be used to improve the shear locking behavior for curved shell configurations by activating the "Reduced integration for membrane shear terms" or "Reduced integration for transverse shear terms" checkbox.  If the thickness of the shell elements are less than 1/10 of the length or height, reduced integration is recommended.  If the thickness of the shell elements is greater than 1/10 of the length or height, full integration is recommended.

 

For Co-rotational shell elements, the "Suppress drilling DOF" checkbox can be activated to help models converge.  This option can not be used if a local coordinate system is being used in the model.


Basic Steps for Using Shell Elements